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We present an algorithm for the explicit numerical calculation of SU(AQ and SL(A^, C) Clebsch-Gordan 
coefficients, based on the Gelfand-Tsetlin pattern calculus. Our algorithm is well suited for numerical imple- 
mentation; we include a computer code in an appendix. Our exposition presumes only familiarity with the 
representation theory of SU(2). 



I. INTRODUCTION 

Clebsch-Gordan coefficients (CGCs) arise when decomposing the tensor product V s ® V of the representation 
spaces of two irreducible representations (irreps) S and S' of some group into a direct sum V Sl ® • • • ® Y Sr of 
irreducible representation spaces. They describe the corresponding basis transformation from a tensor product basis 
{\M <S> M'}} to a basis {|M")} which explicitly accomplishes this decomposition. 

CGCs are familiar to physicists in the context of angular momentum coupling, in which the direct product of two 
irreps of the SU(2) group is decomposed into a direct sum of irreps. SU(3) Clebsch-Gordan coefficients arise, for 
example, in the context of quantum chromodynamics, while SU(AT) CGCs, for general N, appear in the construction 
of unifying theories whose symmetries contain the SU(3) x SU(2) x U(l) standard model as a subgrou 
CGCs are also useful for the numerical treatment of models with SU(iV) symmetry, where they arise when exploiting 
the Wigner-Eckart theorem to simplify the calculation of matrix elements of the Hamiltonian. Such a situation 
arises, for example, in the numerical treatment of SU(A^)-symmetric quantum impurity models using the numerical 
renormalization grourM Such models can be mapped onto SU(iV)-symmetric, half-infinite quantum chains, with 
hopping strengths that decrease exponentially along the chain. The Hamiltonian is diagonalized numerically in an 
iterative fashion, requiring the explicit calculation of matrix elements of the Hamiltonian of subchains of increasing 
length. The efficiency of this process can be increased dramatically by exploiting the Wigner-Eckart theorem, which 
requires knowledge of the relevant Clebsch-Gordan coefficients. (Details of how to implement SU(iV) symmetries 
within the context of the numerical renormalization group will be published elsewhere.) Similarly, tremendous gains 
in efficiency would result from developing SU(iV)-symmetric implementations of the density matrix renormalization 
group for treating generic quantum chain models^', or generalizations of this approach for treating two-dimensional 
tensor network models^. 

For explicit calculations with models having SU(iV) symmetry, explicit tables of SU(AT) Clebsch-Gordan coefficients 
are needed. Their calculation is a problem of applied representation theory of Lie groups that has been solved, in 
principle, long ago^^. For example, for SU(2) RacafP^ has found an explicit formula that gives the CGCs for the 
direct product decomposition of two arbitrary irreps S and S' . For SU(iV), explicit CGC formulas exist for certain 
special cases, e.g. where S' is the defining representation^"!!!. Moreover, symbolic packages such as the program 
"Lie'CEl a i so allow the computation of certain CGCs, but rather have been conceived as a general-purpose software 
for manipulating Lie algebras than a high-speed implementation for calculating CGCs. However, for the general case 
no explicit CGC formulas are known that would constitute a generalization of Racah's results to arbitrary N, S and 
S'. 

The present paper describes a numerical solution to this problem, by presenting an elementary but efficient algo- 
rithm (and a computer implementation thereof) for producing explicit tables of CGCs arising in the direct product 
decomposition of two arbitrary SU(AT) irreps, for arbitrary N. (Since SU(AT) and SL(N, C) have the same CGCs, 
our algorithm also applies to the latter, but for definiteness we shall usually refer only to the former.) Our work is 
addressed at a readership of physicists. Our algorithm uses only elementary facts from SU(N) representation theory, 
which we introduce and summarize as needed, presuming only knowledge of SU(2) representation theory at a level 
conveyed in standard quantum mechanics textbooks. Previous attempts at formulating an algorithm for calculating 
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SU(iV) CGCs are either not sufficiently general for our purpose d 16 * 17 !, or require mathematical methods^ much more 
advanced than ours, far beyond the scope of a standard physics education. 

We begin in Sec.lTTlby formulating the problem to be solved in rather general terms. To set the scene for its solution, 
VII| summarize the various elements of SU(iV) representation theory (without proofs, since this is all 



sections 



III 



to 



textbook material) . First, in Sec. |III| we review the calculation of SU(2) CGCs using a strategy that can readily be 
generalized to the case of SU(iV). Then we proceed to SU(7V) representation theory and review in sections IV to VII 
a scheme, due to Gelfand and Tsetlin (GT^l, for labeling the generators of the corresponding Lie algebra su(N) , 



its irreps and the states in each irrep. The GT-scheme is convenient for our purposes since it yields explicit matrix 
representations for any SU(7V) irrep (Eqs. 
formulate, in sections 



VIII to XII 



28| and (29) below). With these in hand, we are finally in a position to 



our novel algorithm for computing SU(-/V) CGCs: it is simply a suitably generalized 
version of the SU(2) strategy of Sec. |III| 

The main text is supplemented by several technical appendices. App. [X] reviews the relation between the GT- 
patterns used in the text and Young tableaux, with which physicists are perhaps somewhat more familiar. App. IB] 
deals with the Littlewood-Richardson rule for determining which irreps V s " occur in the decomposition V s <g> V . 
App. [C] describes two algorithms, needed for indexing purposes, which map the labels of irreps and of carrier states, 
respectively, onto natural numbers. Finally, App. [dJ which is available in electronic forrrPS, gives the source code for 
our computer implementation, written in CH — h As a service to potential users, we have set up a web site^ containing 
an interactive "CGC-generator" . It allows visitors to perform a number of tasks on input data of their own choice, 
such as finding all irreps S" occuring in the decomposition of S ® S', or finding the complete set of CGCs arising in 
the decomposition of S <£) S' . 



II. STATEMENT OF THE PROBLEM 



To fix notation, let us state the problem we wish to solve for a general matrix Lie group Q . (In subsequent sections, 
we restrict attention to Q — SU(iV) or SL(iV, C).) Let S be an irrep label that distinguishes different irreps of Q of 
SU(JV), and ds the dimension of irrep S. Let V s = span{|Af)} denote the carrier space for S, spanned by d$ carrier 
states \M), where the label M will be understood to specify both the irrep S and a particular state in its carrier 
space. (This will be made explicit in subsequent sections.) Note that, throughout this paper, we adopt the viewpoint 
of quantum mechanics, where we consider only representations on complex vector spaces. Besides, a state is to be 
understood as a one-dimensional subspace, not a vector. However, we pick a representative vector \M) of each such 
subspace and subsequently treat a state as a vector. We assume the inner product of two such normalized vectors 
\M) and \M') to be given by (M\M') — 8m, m' unless noted otherwise. 

The action of a group element g £ G can be represented on V as a linear transformation 



g:\M) 



M' 



(U, 



MM 



\M') 



(1) 



where the Ug are ds X ds dimensional unitary matrices respecting the group structure Ug ± Ug 2 — Ug ig2 - 

Now consider the direct product of two carrier spaces, V <£> V = span{|M ® M')}, of dimension ds • ds>- We are 
interested in its decomposition into a direct sum of carrier spaces V s of irreps S" , 



ee 



rS" 



SS' 



(2) 



Here the integer Ncj s , > 0, called the outer multiplicity of S", specifies the number of times the irrep S" occurs in this 
decomposition, and for a given S" , the outer multiplicity index a = 1, . .. ,N§ S/ distinguishes multiple occurrences of 
S" . Correspondingly, let {\M",a}} be a basis for the direct sum decomposition, i.e. V s ,Q = span{|M", a)}. Carrier 
space dimensions add up according to ds ■ ds> = J2s" ^SS'^S"- 

The decomposition ^ implies that a basis transformation C can be found from the direct product basis to the 
direct sum basis which block-diagonalizes the matrix representations of all group elements (Ref. p. 100): 



C(U, 



9 ® U 9 



5" 



)C f = 



\ 



(3a) 



3 



where each Sj is a shorthand for a certain (S", a). 

Since Q is a matrix Lie group (SU(iV) or SL(iV, C)), it is convenient to work with its associated Lie algebra g (su(7V) 
or si(N, C)). It is obtained by considering the infinitesimal action of Q on V , i.e. by taking derivatives of the group 
at the identity. This derivative acts on the direct product of two group representations according to the product rule, 
so that the basis transformation C could equally be defined by the property that it block-diagonalizes the algebra 
representation: 



c(c/i <g) i s ' +i s ® u$ )c f = 



\ 



\ 



■! 



(3b) 



When projected to the subspace V s ,a (denote the corresponding projector by P s :Q ), the action of the algebra in 
the direct product representation can thus be written as 



c{u%®i s ' +i s '<g)f/!')c ,t p — > 

Concretely, the basis transformation C can be expressed in the form 

\M",a)= £ C%'$\M®M') 



U 



S",a 



(4) 



(5) 



MM' 



where the C M ^, are the Clebsch-Gordan-coefficients of present interest. They are understood to be defined only 

for Ng" s , ^ 0, and express the carrier states of Y s ' a in terms of linear combinations of product basis states from 

V s ® V s . The CGCs encode so-called selection rules, in that C M1 ^ ^ only for a limited number of combinations 
of M, M' and M" . 

Since the CGCs are the entries of the unitary matrix C, they satisfy the following orthonormality conditions: 



° M,M'\ Pm,M' ' ~ °M" ,M"° a >°> ' 



MM' 



E^M" ,ot(f-\M" ,a\* _ c 
^MM'^MM'' ~ 



MM M'M' 



(6a) 
(6b) 



M",a 



Actually, the C M ^ can always be chosen to be real, and we shall do so throughout. 

The goal of the present work is to present (and implement on a computer) an efficient algorithm for Q = S\J(N) 
or SL(7V, C) which, for any specified N and any specified irrep labels S and S' , produces explicit tables of all CGCs 
arising in the direct product decomposition 



III. REVIEW OF SU(2) CLEBSCH-GORDAN COEFFICIENTS 

Before considering the general SU(7V) case, we first review a method for calculating SU(2) CGCs. While there are 
various ways to accomplish this task, the particular approach presented below illustrates the general strategy to be 
used for SU(iV) in later sections. The discussion is structured as follows: First, we recall the Lie algebra associated 
with SU(2), then its irreducible representations, then move on to product representation decompositions, and finally 
set up equations specifying the CGCs. 

The Lie algebra associated with SU(2), denoted by su(2), consists of all real linear combinations of three basis 
elements, J x , J y , and J z , obeying the commutation relation [J x , J y ] = iJ z (plus cyclic permutations of the indices). 
However, it will be more convenient to deal with complex linear combinations of these, which constitute the algebra 
sl(2, C). As a basis for the latter, it is common to choose three elements, J+ = J x + iJ y , J- — J x — iJ y , and J z , 
obeying the following commutation relations: 



[J Z ,J±]=±J±, 
[J+,J-]=2J Z . 



(7a) 
(7b) 
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FIG. 1. SU(2) weight diagram for S — 2. Arrows show the action of J± on the state \S = 2,m = 1). 



Each su(2) irrep, and correspondingly, each SU(2) irrep, can be uniquely (up to an isomorphism) identified by a 
nonnegative half-integer, S — 0, 1/2, 1, . . .. The carrier space V s of such an irrep has an orthonormal basis where the 
states, denoted by \S,m), are labeled by a half-integer, m = S, S — 1, . . . , — S, such that the action of J z and J± is 
given by 

J z \S,m) = m\S, m) , (8a) 
J± \S, m) = v / (S±m + 1)(S T m) \S, mil). (8b) 

The J z -eigenvaluc m will be called the z- weight of the state \S, m) (in anticipation of similar nomenclature to be used 
for SU(iV) below). The action of J± can be visualized in a so-called weight diagram, which represents each carrier 
state \S,m) by a mark on an axis at the corresponding m- value. For example, the carrier space of S = 2 is shown 
in Fig. [I] In anticipation of the generalization to SU(iV), we label basis states from now on by a composite index 
M = (S, 771 ), which includes both the irrep label S and the basis index m. 

Each carrier space V contains a unique (up to normalization) highest-weight state, \H"), defined by the property 
that 

J+ \H) - . (9) 

For su(2), it carries the labels \H) = \S, m — S). 

In the direct product decomposition of two su(2) irreps S and S' , the outer multiplicity Ng s , in the notation of 
Eq. (2| is given by: 

N s" fl for|S-S'|<S"<S + S', 
S ' S ' 10 otherwise. 

oince Nf s , < 1 for su(2), we shall, throughout this section, omit the index a appearing in Eq. ([5]). In particular, 
Eq. (JsJ now takes the form 

\M")= E < M ,\M®M'), (11) 



M.M' 



where the CGCs C$ M , satisfy the selection rule: 



M,M- 

M 



m" ^m + m' => C§ M , = . (12) 

It reflects the fact that \M"), \M) and \M') are eigenstates of jf , jf and Jf , respectively, where the superscripts 
on J z indicate which carrier space the respective operators act on. 

To obtain the CGCs for given S and S' explicitly, we consider each S" for which Ng s , > separately. Let us make 
the following ansatz for the expansion of \H") in terms of product basis states: 

W") = J2 C M " M ,\M®M>), (13) 

where C MM , are the CGCs of \H"), and the sum runs only over values of m and ml that satisfy the selection rule 
Q . Inserting Q into we obtain 



G M,M> (J+ ® I 5 ' + I S ® 4') l M ® M ') = ■ ( 14 ) 



M.M' 
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After evaluating the action of the raising operators on \M <g> M') using Eq. (8b) and requiring the coefficients in front 



of each state \M ® M') to vanish independently, we obtain a homogeneous linear system of equations. We solve for 
Cr, M , and fix a solution by the normalization condition ( 6a ) and by requiring C|> M , to be real and positive for the 



is nonzero. 



largest value of m for which M , 

The CGCs of lower-weight states (i.e. states other than the highest-weight state) are found by noting that 

\M") = \S",m") = N{J-) s "- m " \H") 

=AfJ2 cf; M ,{j s ®i s ' +t 

M,M' 



J: 



S'\S"-m" 



M ® M') 



(15) 



(JV = y/fjS" + m")\/(S" - m")\(2S")\ is a normalization constant 
known from Eq. ( 8b ) . By rewriting it into the form of Eq. ( |11[ ) 



the desired C^m 



The right-hand side of this equation is fully 
can readily be identified. 

For given S", S and S' it is possible to write Eq. (15) as a recursion relation relating CGCs with different m'^l 
Moreover, for su(2), there exists a closed formula for Cjjf Nevertheless, for present purposes, the approach 

presented here is the most convenient as its key steps can readily be generalized to calculate su(N) Clebsch-Gordan 
coefficients. The differences in comparison to su(2) will lie in (i) the more complex structure of raising and lowering 
operators, (ii) the labeling schemes for irreps and states, and (iii) the method for finding the irreps occurring in a 
product representation decomposition, all of which we tackle in the following sections. 



IV. THE LIE ALGEBRA ASSOCIATED WITH SU(iV) 

Instead of working with the group SU(iV) itself, it will be more convenient for our purposes to consider its associated 
Lie algebra, su(iVj23j (ch. 13). The latter consists of all traceless anti-Hermitian n x n matrices, while the ordinary 
commutator serves as its Lie bracket. Most results obtained for representations of Su(N) carry over to SU(iV) one-to- 
one, with the elements of the Lie algebra representing the generators of the Lie group. Notably, the Clebsch-Gordan 
coefficients of their representations are identical. 

We begin by specifying a basis for the su(N) algebra, in order to illustrate its structure. Let E p ' q be the single-entry 
matrices, i.e. = S p>r 6q lS . A possible choice of basis is given by the matrices i(E k ' 1 + E l ' k ) and E k ' 1 — E l,k for 
1 < k < I < N, and i(E 1 ' 1 — E l+1 ' L+1 ) for 1 < I < N — 1. su(N) is spanned by real linear combinations of these 
matrices. Just as for su(2) and sl(2, C), however, it will be convenient to work with a basis for si(N, C). To this end, 
define for 1 < I < N — 1 the complex linear combinations, 

Jjp = *(E l > 1 - E l+1 ' l+1 ), (16a) 

4 ] = E Ll+ \ (16b) 
J W - E l+1 \ (16c) 
which satisfy, for each I, the familiar su(2) commutation relations of Eq. 



± 



= ±4\ (17a) 
= 2J«. (17b) 



The N — 1 matrices form a maximal set of mutually commuting matrices, [Jz"\ ji' '] = (thus, the iJ^ span the 
Cartan subalgebra of su(7V)). Thus, none of the J±* commutes with all elements of this set, or with all other J± ^ 
operators. 

The matrices and J±* are not anti-Hermitian and thus do not belong to su(N), but rather to sl(N, C). However, 
it is sufficient to restrict our attention to Jj/ because, from these, we can recover an anti-Hermitian basis using 

&«=[J^ l \[j!?-* ) ,...[J^ 1) ,jto]]...] for p>q, (18a) 

^ = [J^ ) ,[J^ 1 \...[J^- 2) ,J^- 1) ]]...] for p<q. (18b) 

In other words, once we know representations for all j£ on a given carrier space, the representations of all other 
elements of both the algebras sl(N, C) and su(A^) are also known. For definiteness, we shall refer to su(N) below, 
although the constructions apply equally to sl(N, C). 
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V. LABELING OF IRREPS AND STATES 



The su( N) basis defined in the preceding section has a feature that makes it particularly convenient for our purposes: 
if one also adopts a specific labeling scheme, devised by Gelfand and Tsetlin (GT}^, for labeling su(N) irreps and the 
basis states of their carrier spaces, these basis states are simultaneous eigenstates of all the matrices J z l \ and explicit 
formulas exist for the matrix elements of the J± with respect to these basis states. The next three sections are 
devoted to summarizing the GT labeling scheme without dwelling on its mathematical roots - the mere knowledge of 
its rules is sufficient for our purposes. (The relation of the GT-schemc labeling scheme to a frequently-used alternative 
but equivalent labeling scheme, employing Young diagrams and Young tableaux, is summarized, for convenience, in 
Appendix [A]) 

Up to equivalent representations, each su(N) irrep can be identified uniquely by a sequence of N integers^, 

S = (mi.jv, ■ • - ,mN,N), (19) 

or S = (mfcjv) in short, fulfilling m k ^ > "nik+i.N for 1 < k < N — 1. We shall call such a sequence an irrep weight or 
i-weight, in short. The second index, N, identifies the algebra, su(N); the reasons for displaying this index explicitly 
will become clear below. Two i- weights S and S' for which all components differ only by a /c-independent constant, 
i.e. m' k N = mfcjv + c with c € Z, designate the same su(N) irrep. This fact can be used to bring any i-weight into a 
"normalized" form having ttin.n = 0, which will be assumed below, unless otherwise specified. 

GT exploited the fact that the carrier space of any su(N) irrep splits into disjoint carrier spaces of su(N — 1) irreps 
to devise a labelling scheme with a very convenient property: It yields a remarkably simple rule for enumerating which 
su(N — 1) irreps occur in the decomposition of S = {m k ^), namely all those with i-weights (mi.jv-i, • ■ • , mjv-i,iv-i) 
that satisfy the condition rrik,N > Tn k ,N-i > m/c+i.jv for 1 < k < N — 1. Note that, here, it is crucial not to set 
m N-i,N-i — so that we can distinguish between multiple occurrences of the same su(N — 1) irrep. 

Recursively, the carrier spaces of su(N — 1) irreps give rise to su(N — 2) irreps and so on, down to su(l), the carrier 
spaces of which are one-dimensional. This sequence of decompositions can be exploited to label the basis states \M) 
of a given su(N) irrep 5 = (m,k,N) using so-called Gelfand- Tsetlin patterns (GT-patterns). These are triangular 
arrangements of integers, to be denoted by M = (rrikj), with the structure 



M 



/mi, at m 2 .N ■ ■ ■ m N ^ N \ 

mi jv-l • • • m N-l,N-l 



V mil / 



(20) 



i.e. the first index labels diagonals from left to right, and the second index labels rows from bottom to top. The top 
row contains the i-weight (mk,N) that specifies the irrep, and the entries of lower rows are subject to the so-called 
betweenness condition, 

m kt i > m k ,i-i > m k +i,i (1 < k < I < N). (21) 

The dimension of an irrep S = (to/-,jv) is equal to the number of valid GT-patterns having S as their top row. There 
exists a convenient formula for this number: 

m k .N - m k '.N 



l<k<k'<N 



dim(5)= n I 1 * V-7 J' (22) 



Note that the SU(2) basis state conventionally labeled as \j, m) corresponds to the GT-pattern ( 2j J _ m °), and the 
above formula reduces to dim(j) = 2j + 1. 

To obtain a complete description of SU(7V) irreps, we need to specify how the Lie algebra su(N) acts on states 
labeled by Gclfand-Tsetlin patterns. The following two sections are devoted to this task, section 



section VII dealing with J±K 



VI 



with Jz^ and 



VI. WEIGHTS AND WEIGHT DIAGRAMS 



A very convenient property of the GT-labeling scheme is that every state \M) is a simultaneous eigenstate of all 
JW \M) = \f \M) , (1<1<N-1), (23) 



Jz^ generators 
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with eigenvalues 

Af f = of + afU) (1 < I < N - 1), (24) 

where the row; sum of 1 — X)L=i m fc,i denotes the sum over all entries of row I of GT-pattern M {a^f = by 
convention). We shall call the sequence of TV — 1 jjp eigenvalues the z-weight of the state \M), and denote it by 
W Z (A1) = (A^ , . . . , Ajylj). The z-weight of \M) is a straightforward generalization of the quantum number m in 
quantum angular momentum. 

As will be elaborated below, the notion of weights of states is useful for elucidating the structure of carrier spaces 
of su(iV) irreps, and in particular for visualizing the action of raising and lowering operators. The above way of 
introducing weights is, however, not unique. We shall often find it convenient to employ an alternative definition of 
the weight of states, which has the convenient property that it always yields nonnegative integer elements (in contrast 
to W Z (M)). This alternative weight, to be called pattern weight or p-weight, and denoted by W(M), is defined to be 
a sequence of N integers, W(M) = (wf 1 , . . . , ), where 

wf^ar~ofU (1<1<N) (25) 



is the difference between summing up rows I and Z — 1 of the GT-pattern M. Note that the number of independent 

M 
N ' 



elements of W(M) is the same as that of W Z (M), namely N — 1, since the wf 1 satisfy the relation Yli=i w f I ~ a ¥ 



The two types of weights are directly related to each other: via Eq. (24), we obtain Af 1 — (wf 1 — wM 1 )/2. For 
dcfiniteness, we will mostly refer to p-weights below (noting here that most statements involving p-weights can be 
translated into equivalent statements involving z-weights). 

At this point, the first of several fundamental differences between su(2) and su(N) with N > 3 appears. While for 
su(2), there always exists exactly one state with a given p-weight, this is not the case for su(N) in general; for N > 3, 
several linearly independent states in the carrier space can have the same p-weight. Indeed, two states have the same 
p-weight, W(M) = W(M'), if and only if they have the same set of row sums (af 1 = of 1 for 1 < I < N — 1) (i.e. 
they differ only in the way in which the "weight" of the row sums is distributed among the entries of each row). For 
a given p-weight W, the number of states \M) having the same p-weight, W(M) — W, is called the inner multiplicity 
of that p-weight, to be denoted by I{W). Consequently, p-weights or z- weights are not suited for uniquely labeling 
states of a carrier space (which is why GT-patterns are used for this purpose). 

z-weights nevertheless do provide a convenient way to visualize the carrier space of an su(A^) irrep. To this end, 
consider W Z (M) = (Af 1 , . . . , X^_i) as a vector in (N — l)-dimensional space and, for each state, mark the endpoint 
of its weight vector in an (N — l)-dimensional lattice. The resulting diagram is called a weight diagram. For the su(2) 
irrep j, weight diagrams consist of a coordinate axis with markings at —j, —j + 1, . . . ,j (see Fig. [lj; for su(3), weight 
diagrams are two-dimensional (see Fig. [2]); for N > 4, weight diagrams cannot be readily drawn on paper because the 
corresponding lattices have more than two dimensions. 

Note that, in Fig.jiJ the z-weight W z = (0,0) has inner multiplicity two, since the two states ^% 1 o °j and ^~ i \ °j 
have the same row sums. 



VII. RAISING AND LOWERING OPERATORS 

Weight diagrams are also very convenient for visualizing the action of the raising and lowering operators J± / . The 
action of J± on a given state \M) produces a linear combination of all states of the form \M ± M k ' 1 ) with arbitrary 
k, where this notation implies element-wise addition and subtraction of the single-entry pattern M k ' 1 having 1 at 
position k,l and zeros elsewhere, 



M k,i = 



/0 




V 








/ 



(26) 



(Note that M k ' 1 on its own is not a valid GT-pattern.) Thus the resulting patterns differ from M only in row I. All 
states \M ± M that are generated in this fashion have the same row sums, z-weights, and p-weights (independent 




FIG. 2. Weight diagram of the su(3) irrep (2,1,0). Each dot represents a z-weight; we also indicate the GT-patterns of the 
corresponding states. The double circle around (0, 0) indicates that there are two states with this weight. The solid and dashed 
arrows represent the action of and ji 2 ', respectively. could be represented by arrows pointing in directions opposite 

to those of J^K) Note that both acting on ( 2 ° j and acting on ( 2 1 J produce linear combinations of ( 2 ° J 

and ^ 1 1 J , albeit different ones. (In the literature it is not uncommon to choose a different su(3) basis that renders this 
weight diagram more symmetric.) 



of k), 



M 
1 ' ' 



W z (M±M k ' 1 ) = (A 
W{M±M k > 1 ) = (wf , 



,A 



M \ 

l-2> A l 



M x T 1/2, Xf ± 1, Xflx =f 1/2, A^ 2 ) • 



x M 



(27a) 
(27b) 



unless states with this weight do not exist, in which case the result vanishes. 

The weight-shifting action of lowering operators is illustrated in Fig. [2] for the weight diagram of the su(3) irrep 
S = (2, 1, 0). Since the weight diagram is two-dimensional, there are two lowering operators, J_ and which shift 
in different directions (indicated by solid/dashed lines), (ji produces a shift in the opposite direction of J_ .) Note 
that there are two different "paths" to reach the z-wcight (0,0) from the z-weight (^, |), namely via either J a) ji 2) 

or Since J_ and do not commute, these paths are inequivalent; indeed, they produce two different 

linear combinations of the two states with z- weight (0,0). More generally, the fact that inner multiplicities larger 
than 1 arise for su(N) representations with TV > 2 is a direct consequence of the fact that there are, in general, several 
different ways of reaching one state from another via a chain of raising and lowering operators, and that these ways 

are not equivalent, because j£ and J± ^ do not commute for I V. 

Very conveniently, closed expressions have been found by Gelfand and TsetlkPD for the matrix elements of all raising 
and lowering operators with respect to the basis of GT-patterns. Explicitly, the only nonzero matrix elements of J_ 
are given, for any 1 < k < I < N - 1, by 26 (p. 280): 



(M 



M k ' l \J { l ] \M) 



1+1 1-1 

I! ( m k'.i+i - m k j + k - k' + 1) J] ( m k'.i-i - m k ,i + k- 

k' = l k' = l 

I 

11 i m k',i - mk,i + k - k' + l)(m,k',i - m k j + k - k') 

k' = l 
k'^k 



V 



k') 



(28) 
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These matrix elements are real and nonnegative, and the right-hand side vanishes if M — M k ' 1 is not a valid pattern. 

1 its 

of the preceding formula and replacing \M) by \M + M k ' 1 ): 



As jl is the Hermitian transpose of J_ , we can obtain its nonzero matrix elements by taking the complex conjugate 



(M + M k ' l \J [ l ] \M) 



{ i+i i-i ^ 

il ( m k'.i+i - m k ,i + k-k') J] (mk',1-1 ~ mk,i + k - k' - 1) 

k' = l k'=l 



V 



I 

]1 i m k',i - mk,i + k— k')(mh',i - rrik,i + k- k' - 1) 

k' = l 
k'^k 



(29) 



/ 



These formulae generalize Eq. (8b) to su(N). 

Each irrep has a unique state \H), called its highest-weight state, that is annihilated by all N — 1 raising operators 



JY> Iff) = 



(1 < I < N - 1) 



(30) 



Since \H) is a unique state, the inner multiplicity of its p-weight W(H) is one, and the irrep can be identified by 
specifying W(H). Our labeling scheme indeed exploits this fact: the i-weight of an irrep is equal to the p-weight of its 
highest-weight state \H), i.e. S — W(H). Conveniently, the GT-pattern H = (hk,i) has the highest possible entries 
fulfilling Eq. (21 ), i.e. h k ,i — hk,N io? 1 < k < I < N — 1 (all entries on the k-th diagonal are equal to m k ^)- 



This concludes our exposition of those elements of SU(JV) representation theory in the GT-scheme that are needed in 
this work. In the following sections we discuss the decomposition of direct product representations and the calculation 
of the associated CGCs. The specific details of the strategy described below are, to the best of our knowledge, original. 



VIII. PRODUCT REPRESENTATION DECOMPOSITIONS 

The product of two irreps, say S ® S' , is, in general, reducible to a sum of irreps (Eq. (l2j)). While it is well-known 



for su(2) which irreps occur in such a decomposition (see Eq. (10)), the corresponding result for su(N) relies on a 
relatively simple but hard to prove method based on the Littlewood- Richardson rule^. This method involves writing 
down all possible GT-patterns for the irrep S and using each of these to construct, starting from 5", a new irrep S" . 
As the outcome of this method is the same when interchanging S and 5", it is preferable to take the irrep with the 
smaller dimension of the two as S. 

For given irreps S = (m^jv) and S' — (m' k N ), and a particular GT-pattern M = (wikj) associated with S, let us 
introduce some auxiliary notation. For I = 1, . . . , N and k = 1, . . . , I, we set 6^; = m,k,l ~ m k,i-i (where rrik i = if 
k > I, for case of notation) and B^ i = m\ N + Ylk'=i (note that here, m' k l carries a prime, while b k i does not). 
Then, the irrep S" — (m'l N ) = (Bk.k) occurs in the decomposition of 5* ® S' if and only if 

Bk-i,i > Bfe-1,2 > • • • > > B kt i > B kM1 >■■■> B k , N for all 1 < k < I < N . (31) 



(We emphasize that this condition must hold for each value of k and I.) By checking whether (31) holds for all 
GT-patterns associated with S, all 5" in the decomposition of S (8> S' can be identified. 



There exists a more efficient way to validate Eq. (31 ) than to check each value of k and I independently. For a given 



GT-pattern M = (mk,i) associated with S, proceed as follows: 

1. Initialize (ti, . . . , tjy) — (m! x N , . . . , m' N N ) by the i-weight of S' . 

2. Step through the pattern M along the diagonals from top to bottom and from left to right, i.e. in the order 
mi,iv, m^AT-i, • • ., mi,i, m 2 ,N, ma.jv-i, ■ ■ ■, w 2 ,2, • ■ •, rn N . N . 

3. At each position, say n%k,i, replace ti by t\ + b k j. 

4. If I > 1, check whether f;_i > ti- If this condition is violated, discard this GT-pattern, construct the next one, 
and commence again from step 1. 

5. If we reach the end of the pattern M, the current value of (ti, . . . , %) specifies the weight of an irrep S" that 
occurs in the decomposition of S <8> S' . 
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For N > 2, this procedure in general can produce several occurrences of the same irrep S" . The number of such 
occurrences, denoted by Ng S , in Eq. ([2]), is the outer multiplicity of S". (For SU(2), the outer multiplicity is either 
or 1.) 

Let us illustrate this procudure by an example (individual steps are shown in Table |l| : 



(2, 1, 0) ® (2, 1, 0) = (4, 2, 0) 8 (3, 3, 0) © (4, 1, 1) © (3, 2, 1) © (3, 2, 1) © (2, 2, 2) 
= (4, 2, 0) © (3, 3, 0) © (3, 0, 0) © (2, 1, 0) © (2, 1, 0) © (0, 0, 0) 



(32a) 
(32b) 



(For the second line, we adopted "normalized" i-weights with m N N — 0.) To check that the dimensions are correct, 



use Eq.(22 1 to verify the dimensions of the irreps in this equation are 8x8 = 27 +10 + 10 + 8 + 8 + 1, respectively. 



Note that the irrep (2, 1,0) occurs twice in the decomposition, in other words, its outer multiplicity is 2. 



IX. SELECTION RULE FOR SU(JV) CLEBSCH-GORDAN COEFFICIENTS 



The fact that all states labeled by GT-patterns are eigenstates of operators implies a selection rule for SU(iV) 
CGCs. Explicitly, let us consider a state \M") occurring in a decomposition of a product representation. On the one 
hand, we have 



J<P\M",a) = Xf"' a \M",a), 



(33a) 



and on the other hand, by Eqs. Q and (JSJ, 

J« \M", a)=J2 Cm'm'(4 1) ' S ® I 5 ' + I S ® 4 l) ' S ') \M ®M')=J2 C%" t £{\ M + xf) \M © M') . (33b) 



MM' 



MM' 



These equations can only be fulfilled if vanishes whenever Xf 1 ,a + 1 Xf 1 + Xf ! for any I. Defining an element- 

wise addition on weights, we write, in short: 



W Z (M") ? W Z (M) + W Z {M') => C%'£, = 



(34) 



This equation (or a transcription thereof involving p-weights) represents the generalization of Eq. ( 12 ) to su(iV) 



CLEBSCH-GORDAN COEFFICIENTS OF HIGHEST-WEIGHT STATES 



After determining which kinds of irreps S" appear in the decomposition of a product representation, we are ready 
to construct their Clebsch-Gordan coefficients. For each S" , we start by finding the CGCs of its highest-weight state, 
\H" , a), as defined in Eq. (30). The index a — 1, . . . , Ng s , distinguishes between the instances of irreps with outer 
multiplicity. Nevertheless, we determine the CGCs of \H",a) with given S" for all values of a in a single run. 

For this purpose, we make an ansatz of the form §5§ for the highest-weight state (compare Eq. (13)), 



\H",a)=Y;C M "£,\M®M>) 

MM' 

W{M) + W(M') = W{H" xt) 



(35) 



with CGCs C^j f?,, where the sum is restricted to those combinations of states \M © M') that respect the selection 



rule (34). Now insert Eq. (351 into Eq.(30) to obtain (compare Eq. (14)), 



E r H", a( T (l),S 

MM' 

W(M)+W(M')=W(H" ,a) 



< I s ' + I s © jj 0,5 ') \M © M') = 0, 



(1 < I < N - 1). 



(36) 



After evaluating the action of the raising operators on the product basis states via Eq. ( 29 ) , we obtain a homogeneous 



linear system of equations in the CGCs C^m'- It has N§ S i linearly independent solutions, one for each value of 
a. Thus, an outer multiplicity larger than 1 leads to an ambiguity among the CGCs of the highest-weight states 
of all irreps of the same kind 5"': a unitary transformation \H,a) — > Yla' U a ,a> \H,a'} among the highest-weight 
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final irrep 


(4,2,0) 

(3,3,0) 
(4,1,1) = (3,0,0) 
(3,2,1) = (2,1,0) 
discarded 

(3.2.1) = (2,1,0) 
discarded 

(2.2.2) = (0,0,0) 


previous + 
(0,0,63,3) 


(4,2,0) 
(3,3,0) 

(4,1,1) 
(3,2,1) 

(3,2,1) 

(2,2,2) 


previous + 
(0,62,2,0) 


(4,2,0) 
(3, 3,0) 
(4,1,1) 
(3,2,1) 

(3,2,1) 

(2,2,2) 


previous + 
(0,0,62,3) 


(4,1,0) 
(3,2,0) 

(4,1,1) 
(3,2,1) 

(3,1,1) 
(3,1,2) 
(2,2,2) 


previous + 
(61,1. 0,0) 


(4,1,0) 
(3,2,0) 
(4,1,0) 
(3,2,0) 

(3,1,1) 
(3,1,1) 
(2,2,1) 


previous + 
(0,6i,2,0) 


(2,1,0) 
(2,2,0) 
(2,1,0) 
(2,2,0) 
(2,3,0) 
(2,1,1) 
(2,1,1) 
(2,2,1) 


'> <=> 

£ 

ft"—' 


s s 1—1 1—1 i—l 

1— (l— <!— <1— 1 1— 1 T 1 T 1 T 1 

of o-f oT cf of of of c£ 


11 3 

> ^ CM 

It 6 

CO 


(O'T'S) 
(O'T'S) 
(O'T'S) 
(O'T'S) 
(O'T'S) 
(O'T'S) 
(O'T'S) 
(O'T'S) 


■■— 1 CO 

a M - 

O CO CM r— 1 

Q. CM rH 

*° * 

S— 1 co 1— ' 
Sh "-o 

8 * 


00000000 

^I^IOOO^OO 

O CNO 1 1 — 1 r-f^H OO 1 — 1 1 — I O 

o^ho^hcmoo^i 

OOOOO^h^h^h 


CO 

CD ^~ 

ft E * 

03 21 


00000000 

,-1^1000^00 

r-f CN^-f ^1 CNl^ ^H^H O^H r- 1 ,H O 

CMCMCMCMCMCNtNCM 



TABLE I. Application of the Littlewood-Richardson rule according to the steps of Section |VIH| for the decomposition of 
S(g> S' = (2,1,0) <g> (2,1,0). 
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states will produce different, but equally acceptable highest-weig ht CGCs The full set of CGCs of the irreps 

S" will change accordingly, too. For some applications, there is no need to uniquely resolve this ambiguity. For 
applications where it must be resolved, we will adopt the following convention, suggested by G. ZarancPM Write down 



the independent solutions in the form of a matrix with elements C 



H",, 

MM' ) 



where a = 1, 



1 N§ s , serves as row index 



and (M, M') = 1, . . . , 1(H) as composite column index (where 1(H) is the inner multiplicity of W(H) in the product 
representation). Then use Gaussian elimination to bring this matrix into a normal form, namely the reduced row 
echelon form, 



C 



H" ,a 
M,M' 










Vo 



+ 











+ 









* 
* 
* 



* 
+ * 



(37) 



where + and * denote positive and arbitrary matrix elements, respectively. This normal form is the same for all 
equivalent matrices. To obtain orthonormal highest-weight states, we then do a Gram-Schmidt orthonormalization of 
the rows of the resulting matrix from top to bottom. This procedure uniquely specifies the CGCs for the highest-weight 
states. 

As an aside, we note that the abovementioned ambiguity does not arise for the case of S' = (1,0,... ,0) (the 
defining representation of su(N)) and arbitrary S, since then all outer multiplicites are either zero or one, i.e. then 
N§ s , = or 1. (However, there would still be a sign ambiguity for the CGCs, and the above procedure constitutes 
one way of fixing it.) We note that for this case, explicit formulas for SU(iV) CGCs can be founcW 



XI. CLEBSCH-GORDAN COEFFICIENTS OF LOWER-WEIGHT STATES 



Let us now turn to the CGCs of states of S" other than its highest-weight state. These are obtained by acting on 
both sides of Eq. (35) with lowering operators, using Eq. (28) for the matrix representations of J_ for the carrier 
space V s ' a on the left-hand side, and for the direct product carrier space Y s ® V s on the right-hand side. However, 



according to Eq. (28), the action of J_ in general produces not a unique basis state, but a linear combination of 



basis states of V s ' a . We shall therefore calculate, in parallel, the CGCs of all basis states with a given a and given 
p-weight W = (wi), i.e. of all \M", a) having W(M") = W. 

To this end, assume that we have already determined all "parent states" of the desired p-weight W within V s ' a . 
By parent states we mean those which, when acted upon by a single J_ , yield (linear combinations of) states of 
weight W. For a given J^p (with 1 < I < N — 1), the relevant parent states have p-weight (wi, . . . , Wi + 1, wi + \ — 
1, wi + 2, ■ ■ ■ , wn) and consist of all states of the form \M" + M k,l ,a) with W(M") = W and 1 < k < I, for which 
M" + M k ' 1 is a valid GT-pattern. Each parent state can be expressed as 



\M" 



M k 'Ka) 



M,M' 



M" +M 
M,M' 



\M <g> M') 



where the CGCs are, by assumption, already known. Now, the action of J_ on any parent state can be written as a 
linear combination of all states \M"',a) with W(M"') = W, 



(38) 



jM // + M k,l ja) = J- \M"', a) , (39) 

M<" 

where the coefficients bji" k \ are determined by the matrix representation of J® within V s ,Q , as given by Eq. (28). 



Combining Eqs. (38) and (39) and using the direct product representation of on V 5 <£> V s , we obtain a linear 
system of equations of the 

(40) 



orm (compare Eq. ( |15[ )): 

]T C'" M \M"', a) = J2 CM,M> Mh ' l ' a (J^' S ® 1 S ' + 1 S ® J^ 3 ') \M ® M') 



M" 



M,M' 
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Each combination of indices M", k, and I specifies a separate equation, where M" runs over all GT-patterns such 
that W(M") = W, I runs from 1 to N - 1, and k runs from 1 to I, provided that M" + M k ' 1 is a valid GT-pattern. 
Actually, only I(W) of these equations are linearly independent; as we do not know in advance which ones these are, 
we include them all, i.e. the system of equations (401 is, in general, overdetermined. Since the action of the J_ s on 



the right-hand side is known from Eq. ( 28 ) , the sought-after CGCs A ; f " can now be readily obtained by inverting 



the matrix of the coeffcients b ^ ,, k ; in order to bring Eq. ( 40 ) into the familiar form of Eq. ^ . 

XII. ALGORITHM FOR COMPUTER IMPLEMENTATION 

Having gathered in the preceding sections all necessary ingredients, we are now ready to formulate the sought-after 
algorithm for calculating SU(7V) CGCs. Given two SU(iV) irreps S and S' , perform the following steps: 



1. Find the irreps S" appearing in the decomposition of S <g> S", as described in Sec. VIII 



2. For each irrep S", find the Clebsch-Gordan coefficients of the N§ s , highest- weight states \H",a). Resolve outer 
multiplicity ambiguities, as described in Sec. |Xj 

3. From each highest-weight state \H",a), construct the lower-weight states by repeated application of J_ oper- 



ators, treating each weight of S" separately, as described in Sec. XI 



An explicit computer implementation of this strategy is presented in App. [TJ] To check that our algorithm works 
correctly, we have verified that it satisfies the following consistency checks: 

• For SU(2) and SU(3), the results coincide with known formulas and tables, up to sign conventions. 



• The selection rule (34 1 is fulfilled. 

• The matrix C of Clebsch-Gordan coefficients (see Sec. [Tl| is unitary. 

• The matrix C block-diagonalizes the representation matrices (Eqs. Q). 

The speed of the algorithm depends polynomially on the dimensions of the irreps S and S". On a modern computer 
(2 GHz CPU clock speed), smaller su(3) cases (e.g. dim S = 6, dim S' = 15) run instantly, while medium-sized su(5) 
cases (e.g. dim S = 35, dim 5' = 224) take a few minutes, and larger su(5) cases (e.g. dim S — 280, dim S' = 420) 
require several hours computing time. 

As an outlook, we note that it should be possible to greatly speed up our algorithm by exploiting the fact that 
the weight diagrams are symmetric under the Weyl group, which in this context can be thought of as the group of 
all permutations of the elements of the p- weights, (wf 1 , . . . , ) — > (w*[^, . . . , w^ N ^). Exploiting this symmetry is 
a nontrivial task, since the Gelfand-Tsetlin basis is not stable unter the operation of the Weyl group. Nevertheless, 
we expect that it should be possible to do within the general framework of our algorithm, by adopting a suitably 
modified state labeling scheme that exploits the Weyl symmetry Work along these lines is currently in progress. 

Acknowledgements: This work was inspired by the progress made by G. Zarand, C. P. Moca and collaborators^ 
in devising a flexible code for the numerical renormalization group, capable of exploiting non-Abelian symmetries. We 
acknowledge stimulating and helpful discussions with G. Buchalla, R. Helling, M. Kieburg, P. Littelmann, C. P. Moca, 
A. Weichselbaum, and G. Zarand, and financial support from SFB-TR12. 

Appendix A: Correspondence between Gelfand-Tsetlin patterns and Young tableaux 

There exists a one-to-one correspondence between i-weights and Young diagrams, and between GT-patterns and 
semi-standard Young tableaux. Thus, our algorithm could equally well have been formulated in terms of Young 
diagrams and Young tableaux. Since the latter are easy to visualize and are perhaps more widely known in the 
physics community than the GT-scheme, this appendix summarizes the relation between the two schemes. Our 
reason for preferring GT-patterns to Young tableaux lies in the complexity of the computer implementation: GT- 



patterns can be stored in a simpler data structure and allow for a simpler evaluation of the matrix elements ( 28 ) 



and (291 



Note that Young tableaux can also be used to label bases that differ from the GT basis used in this work, notably 
the one constructed via Young symmetrizers^ (ch. 7). Thus, the correspondence between GT-patterns and Young 
tableaux set forth below is purely of combinatorial nature. 
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(b) 



1 


1 


1 


1 


1 


2 
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2 
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1 


3| 


2 


2 


2 
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2 
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2 




3 
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3 




3 





TABLE II. | (a) | Examples of Yo ung diagrams of su(3) irreps. Since columns with 3 boxes can be deleted, the last example is 
effectively equal to the first one. (b) Set of all of valid su(3) Young tableaux of shape pp. 



Row 1 



(a) 



V 2 / 
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Row 2 



3 2 
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Row 3 



Row 4 



3 2 1 
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3 2 
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Diagonal 2 
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2 
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Diagonal 3 Diagonal 4 
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V 2 / 
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3 2 1 
3 2 

2 / 
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TABLE III. Conversion of a GT-pattern to a Young tableau, stepping (a) along rows, and (b) along diagonals 



1. Definition of Young diagrams and Young tableaux 

A Young diagram is an arrangement of boxes in rows and columns conforming to the following rules: (YD.l) there 
is a single, contiguous cluster of boxes; (YD. 2) the left borders of all rows are aligned; and (YD. 3) each row is not 
longer than the one above. 

Note that the empty Young diagram consisting of no boxes is a valid Young diagram. For the purpose of describing 
an su(N) irrep, we additionally require that (YD. 4) there are at most N rows; and (YD. 5) columns with TV boxes are 
dropped, i.e. diagrams which differ only by such columns are identified with each other. 

Every Young diagram D satisfying rules (YD.l) to (YD. 5) uniquely labels an su(iV) irrep (or sl(N,C) irrep), i.e. 
the label S used in the main text can be associated with a Young diagram D. Some su(3) examples are shown in 



Table Ha A further example is given by the Young diagrams specifying su(2) irreps: The irrep S — j (describing 
total angular momentum j) corresponds to a Young diagram with 2j boxes in a single row. 

A (semi-standard) Young tableau is a Young diagram, of which the boxes are filled according to the following rules: 
(YT.l) Each box contains a single integer between 1 and N, inclusive; (YT.2) the numbers in each row of boxes 
weakly increase from left to right (i.e. each number is equal to or larger than the one to its left); and (YT.3) the 
numbers in each column strictly increase from top to bottom (i.e. each number is strictly larger than the one above 
it). 

The basis states of an su(N) representation identified by a given Young diagram D can be uniquely labeled by 
the set of all associated valid semi-standard Young tableaux (satisfying rules YT.l to YT.3), i.e. the label M used 
in the main text can be associated with a valid Young tableau T. We shall denote the corresponding state by \T). 



For example, all eight Young tableaux for the diagram CD with respect to su(3) are shown in Table lib As another 
example, let us give the correspondence between states \S, m) of ansu(2) irrep and Young tableaux: \S, m) corresponds 
to a Young tableau with 25* boxes in a single row, containing 1 in the leftmost S + m boxes and 2 in the remaining 
S — m boxes. 

The dimension of a carrier space labeled by a Young diagram is given by the number of valid Young tableaux with 
the same shape as the Young diagram. 



2. Translating GT-patterns to Young tableaux 

Each GT-pattern M = (mkj) uniquely specifies a corresponding Young tableau (p. 526 of Ref. |3D|, which can be 
constructed as follows. Start with an empty Young tableau (no boxes at all), and step through the entries of the 
pattern using either of the following two stepping orders, illustrated in Table [TTla| and pTb"l respectively: 

(a) Proceed from the bottom to top, one row at a time (increasing / from 1 to n), and within each row from left to 
right (increasing fc from 1 to I); or 

(b) Proceed from left to right, one diagonal at a time (increasing k from 1 to n), and within each diagonal from 
bottom to top (increasing I from fc to n). 
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For each step to a new entry in the pattern, say rrik.i located in diagonal k and row I, extend the length of the A:-th 
tableau row to a total of nik.i boxes, by adding to its right boxes containing the number I. 

According to the above procedure, the topmost row of the GT-pattern specifies the number of boxes in the rows of 
the corresponding Young diagram: for the latter, row k of the latter contains mk t N boxes. In this way, the information 
specifying the irrep S, which for a GT-pattern resides in its topmost row, specifies the shape of the corresponding 
Young diagram. Moreover, the number of Z-boxes (i.e. boxes containing the number I) in tableau row k, say dk,i, is 
given by 

dk,i = m,k,i - mk,i-i , (where mu,i = if k > I). (Al) 

Since both stepping orders ensure that pattern entries in the same diagonal k are visited in order of increasing I, 
they yield the same final Young tableau. Order (b) has the feature that an entire tableaux row is completed before the 
next row is begun. As a result, (b) is more convenient for transcribing the Littlewood- Richardson rule for decomposing 
a product representation from the language of Young tableaux to that of GT-patterns. 

The converse process of transcribing a Young tableau to a GT-pattern can be achieved by using the tableau's A:-th 
row, read from left to right, to fill in the pattern's fc-th diagonal, from bottom to top, in such way as to respect the 
above rules. 



3. Remarks about Young tableaux 



In order to aid our intuition for the su(N) representation theory presented in the main text, this section restates 
some of the properties discussed there in terms of Young tableaux. 



The p-weight W(M) of a GT-pattern M, as introduced in Sec. VI has an illustrative interpretation; wf 1 is the 



number of Z-boxes (i.e. boxes containing I) in the tableau corresponding to M. Thus, for the highest-weight Young 
tableau (the GT-pattern of the corresponding state \H) is specified at the end of Sec. VII), row I from the top 

Furthermore, if the states |T) and \T') have the same p-weight, 



,M 



mi,N), e.; 



2 2 

3 3 



contains only Z-boxes (i.e. wf 

the tableaux T and T" contain the same set of entries (i.e. the same number of Z-boxes), but arranged in different 
ways. For example, for su(3) [S3 and S3 have the same p-weight W — (1, 1, 1). 



The action of the raising and lowering operators J± on p- weights is given by Eq. (|27[) . The corresponding action of 



Jj/ on a state labeled by a Young tableau T produces a linear combination of states labeled by tableaux containing 

one more l-box and one less (I + l)-box. Analogously, J_ has the reverse effect on Young tableaux: it produces a 
linear combination of tableaux containing one less Z-box and one more (/ + l)-box. 



Appendix B: Derivation of our formulation of the Littlewood- Richardson rule 

Our formulation of the Littlewood-Richardson rule in Section VIII is based on a version by van LeeuwerP^, formu- 



lated in terms of Young tableaux, which we outline here. We then rephrase this in the language of Gelfand-Tsetlin 



patterns to derive the method presented in Sec. VIII in particular Eq. (31 1. 

Given two Young diagrams D and D' , write down all possible semistandard Young tableaux for D, and for each 
such tableau (to be called the current tableau below) , construct a corresponding Young diagram (to be called the trial 
diagram below) in the following manner: 

1 . Start the trial diagram as a fresh copy of D' . 

2. Step through the boxes of the current tableau from right to left, from top to bottom. 

3. If the box encountered at a given step is an Z-box, add a box at the right end of row I of the trial diagram. 

4. If this produces a trial diagram that is no longer a valid Young diagram (having a row longer than the one 
above), discard it and start anew with the next tableau. 

5. If, however, a valid Young diagram is constructed during each step, the final Young diagram obtained after the 
last step represents an irrep occurring in the decomposition of D <E) D' . 

Let us now translate the above steps into the GT-scheme, thus deriving the rules set forth in Section |VIII| There, 
we assume two i-weights S and S" to be given instead of two Young diagrams. Naturally, taking a fresh copy of D' 
corresponds to initializing (ti, . . . , tjv) = (™'i at, ■ ■ ■ , m'^ at), and stepping through the current tableau in the reading 
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P(S) 


s = 


(mi,4, m 2 ,4, ms.4, m4, 4 ) 







(0,0,0,0) 


1 




(1,0,0,0) 


2 




(1,1,0,0) 


3 




(1,1,1,0) 


4 




(2,0,0,0) 



P(S) 


5 = 


(mi,4,m2,4,m 3 ,4,m4,4) 


5 




(2,1,0,0) 


6 




(2,1,1,0) 


7 




(2,2,0,0) 


8 




(2,2,1,0) 


9 




(2,2,2,0) 



TABLE IV. The first few i-weights of su(4) (excluding weights with 7714,4 7^ 0), arranged in increasing order. 




1 2 3 



di d.2 ■ ■ ■ (ijv-i ^-n 

* ft ft ft t 

m k .N rrit.N -1 + N 



(b) 



FIG. 3. Enumeration scheme of i-weights. | (a) [ Illustration of the combinatorics underlying Pk(S). |(b)| Striking out items such 
that each fhj, ^ takes on the largest possible value. 



order of step 2 corresponds to stepping through the GT-pattern M associated with S along the diagonals from top 
to bottom and from left to right. (This follows from the rules for translating GT-patterns to Young tableaux given 
in Sec. |A 2| recall that the k-th diagonal of a GT-pattern specifies the content of the fc-th row of the corresponding 
Young tableau.) 

Instead of processing one box of the current tableau at a time, we treat all identical boxes of a given row at once 
when stepping through the corresponding GT-pattern. Recalling that bf. 1 of Eq. (Al) gives the number of Z-boxes in 

row k of the current tableau, it follows that B^ i = ml, N + Y^ik'=i then gives the number of boxes in row I of the 
trial diagram after having processed all boxes of type I in row k of the current tableau. 

The condition (31 ), which must be fulfilled for all 1 < k < I < N, finally assures that the trial diagram is a valid 
Young diagram after each step. 



Appendix C: Identifying irreps and states by a single integer 

For numerical codes dealing with i-weights, it is useful to identify each i-weight by a unique number. To this end, 
we need a one-to-one mapping between the set of all su(A^) i-weights (for given N) and the set of nonnegative integers. 
We shall construct such a mapping by devising an ordering rule for i-weights, using this rule to arrange all possible 
diagrams in a list of increasing order, and labeling each i-weight by its position in this list. 

Similarly, we would like to map GT-patterns to matrix indices, so we also need a one-to-one mapping between the 
set of all GT-patterns belonging to a given irrep and the integers from 1 to the dimension of that irrep. Therefore, 
we also define an order on GT-patterns of a given irrep and proceed analogously. 



1. Identifying i-weights with a single number 

We adopt throughout the convention for an i-weight S = (mfc,jv) that tojv.at = (Sec. |v|). 

For i-weights we choose the following ordering rule: the "smaller" of two i-weights is taken to be the one with the 
smaller first element; in case of a tie, compare the second element, and so on. Formally, given two i-weights S and 
S", we assign the order 

S < S' if and only if , for the smallest index (say k) for which m,k,N 7^ m'k.Ni we have mk,N < Tn'k.N ■ (CI) 



Table IV shows the first few i-weights of SU(4), arranged in increasing order. 
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Using this ordering rule, all possible SU.(N) i-weights can be arranged in a list of increasing order and uniquely 
labeled by a nonnegative integer, say P(S), giving its position in this list, 



P(S) = #{S'\S' < S}. 



(C2) 



To determine P(S) for a given i-weight S, we simply count the number of smaller weights S': this number is given 
by the number (say Pi(S)) of all weights S' with m' x N < mi,jv 5 plus the number of all S' with m[ N — mi,jv but 



/??' 



2,N 



< m 2,N ( sa y Pi(S))] etc. Thus, 



N-l 



P(S) = £ Pk(S) , 



(C3) 



fc=i 



where P k (S) is the number of weights S' whose first k — 1 entries are the same as those of S (m' k , N — m k i ^ for 
all k' < k), while the fc-th entry is arbitrary but smaller than that of S (m' k N < to^jv), and the remaining entries 
arbitrary (but subject to 5" being a valid i-weight, with m! N N — 0). The nontrivial "free" (though constrained) 



entries of S', namely (m' h 



k,N' ' 



whose entries rhz ft 



k-l + k,N 



k+1 N , . . . , m' N _ 1 N ), can be viewed as an i-weight S 
(for 1 < k < N) satisfy 



1 > rh x fj > rh 2 „ > 



(m~ k of length N = N - k, 



(C4) 



Pk(S) thus is the number of allowed weights S that satisfy (C4| 



To calculate Pk(S), we note that it is equal to the number of ways to draw or "strike out", from the set of integers 
{1, ... ,nik,N — 1 + N}, an ordered subset {d k } of N integers, d\, < d% < ■ ■ ■ < dft (see Fig. 3a), since there is a 



one- to-one correspondence between the set of all possible such strike-outs and the set of all i-weights S satisfying 



(C4): for a given struck-out set {d k }, with 1 < k < N, set rh k ft equal to the number of non-struck-out integers 



smaller than dft +l _ k (i.e. fh k ft = dft 



JV+l-fc 



(JV + 1 



k)). For example, the weight S that is largest (w.r.t. to the 
ordering rule ((Clj)), namely having all elements equal to m k ^ N — 1, is obtained by choosing the struck-out integers d k 



to be as large as possible (see Fig. 3b I. Thus, we have 



P k (S) 



m k ,N -1 + N' 
N 



N — k + rrik.N — 1 
N-k 



and, consequently, 



P(S) 



E A 1 fN - k + mk,N - 1 
V N-k 

fe=i v 



(C5) 



(C6) 



2. Mapping of Gelfand-Tsetlin patterns to matrix indices 

In analogy to the ordering we have defined on i-weights, we introduce an ordering on the set of Gelfand-Tsetlin 
patterns of a given irrep (i.e. given top row of the pattern). Let M = (rrik,i) and M' = (m' k ; ) (where 1 < k < I < N) 



denote two patterns with m^ jv = m 'k n f° r ^ = 1, - - ■ , ^V". We define a row- by-row ordering of indices (see Table Va|, 
increasing from left to right within a row, and from top row to bottom row, i.e. (k, I) < (k' , V) if I = I 1 and k < k' , or 
if / > I'. We then define M' < M if and only if for the smallest index for which m' k , ^ rrik,i, we have m' k l < rrik,i- 
An example of this ordering is given in Table |Vb| 

We map each Gelfand-Tsetlin pattern M to a nonnegative integer Q(M) by counting the number of smaller Gelfand- 
Tsetlin patterns, i.e. 

Q(M) = #{M'|M' < M} . (C7) 

This number can be determined by generating the pattern (say M({rhk,i})) located directly preceding M in the 
ordered list of patterns, then the pattern preceding M, and so on, until we arrive at the beginning of this list. To 
construct the predecessor of the pattern M, we start by finding the largest index (fc, I) whose entry m k j can be 



decreased without violating the betweenness condition (21 1, rewritten here as 



m k ,i+i > m k .i > m k+u+1 (l<k<l + l<N), (C8) 
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12 3 
4 5 



Q(M) 



( ' 1 1 1 




( 2 1 1 




(2 1 1 




( 2 1 q\ 




^2 1 o\ 




(2 1 q\ 




(2 1 <0 




^2 1 o\ 




< 


1 


< 


1 1 


« 


2 


< 


2 


< 


2 


« 


2 1 


< 


2 1 


V o J 




V i J 




V 1 ) 




V j 




V 1 J 




v . ; 




v x ; 




V = J 


1 




2 




3 




4 




5 




6 




7 




8 



TABLE V. (a) Illustrating the row-by-row rule cho sen in App. C 2 to define an ordering scheme for the indices of GT-patterns: 
(k, 1) < (k', V) if I = V and k < k' , or if I > I', [(b)] Ordering of all GT-patterns belonging to the SU(3) irrep (2, 1, 0), together 
with the corresponding pattern indices Q(M). 



with respect to smaller indices while disregarding it with respect to larger indices (i.e. without violating the second 
inequality, but disregarding the first). Thus, {k, I) is the index for which mk,i — ^fc+i.z+i for all (k,l) > {k,l) but 
m k 1 ^ m k+i i+v We- then decrease j by one and reset the entries of all larger indices to the maximal values that 
satisfy the new betweenness condition. Concretely: 

mf-,i for (k, I) < (k, I) (keep entries with smaller indices unchanged) 
fhk,l — \ rrik^i — 1 for (k, I) = (k, I) (decrease by 1 the entry with largest index for which this is possible) (C9) 
rhk,i+i for (k, I) > (k,l) (give entries with larger indices their largest possible value). 

The number Q(M) is, of course, the number of times we can repeat the process of constructing a preceding pat- 
tern. This procedure maps the lowest-weight and highest-weight states of an irrep S to the numbers 1 and dim(S'), 
respectively. 



Appendix D: Source code 



Below, we provide a C++ implementation of our algorithm, consisting of four fundamental classes: 1. weight is 
a data structure for irrep and pattern weights, 2. pattern stores GT patterns, 3. decomposition implements the 
Littlewood-Richardson rule, and 4. coefficients computes and stores the actual CGCs. The end of the source code 
contains examples of typical applications. For example, to calculate CGCs, perfom the following steps: 

1. Create two objects clebsch: : weight S and clebsch: : weight Sprime, representing the irreps S and S' . 

2. Create the object decomp as clebsch: : decomposition decomp(S, Sprime); this generates the irreps S" that 
occur in the decomposition of S (& S' according to the Littlewood-Richardson rule. (Its output can be read out, 
if desired, as follows: Read out the total number of irreps S" by calling decomp . size () . Read out the j-th one 
of these (with 1 < j < decomp. size ()) by creating an object clebsch: : weight Sdoubleprime (decomp ( j ) ) . 
Read out its outer multiplicity by calling decomp .multiplicity (Sdoubleprime) .) 

3. Pick one of these irreps Sdoubleprime and create the object C as clebsch: : coefficients C (Sdoubleprime , 
S, Sprime); this generates all CGCs C^m? needed for constructing the irrep S" , with multiplicity index a, 
from S and S'. 

4. The Clebsch- Gordan coefficient C^j^f is then read out as C( alpha, Qdoubleprime , Q, Qprime) , where alpha 
indexes the outer multiplicity of S" , and Q, Qprime, and Qdoubleprime are the pattern indices of M , M' and 
M". 

Other common applications involve the translation between an i- weight S and its index P(S), or between a GT-pattern 
M and its index Q(M). To obtain the i-weight index P(S) from the object clebsch: :weight S, call S . indexO , and 
to obtain the pattern index Q(M) from the object clebsch: :pattern M, call M.indexO. Conversely, to construct 
an i-weight S = (mk,N) from a given irrep index P, create the object clebsch: : weight S(N,P), and read out the 
elements mk,N as S(k). Similarly, to construct a pattern M — (rrik,i) in irrep S from a given pattern index Q, create 
the object clebsch: : pattern M(S,Q), and read out the elements rrik i as M(k,l). Finally, to find the dimension ds 
of the irrep S, create the object clebsch: :weight S and call S . dimensionO . 

All of these applications are elaborated in the sample routine main at the end of the source code (starting around 
line 1000). They are also implemented in the interactive CGC-generator available at |http : / /h omepages . physik . | 
uni-muenchen . de/~vondelf t/Papers/ ClebschGordan/ 

To locate the implementation of key equations of the algorithm in the source code, search for the following equation 
numbers: i-weights S: Eq. (19); GT-patterns M: Eq. (20); irrep dimension dim(S'): Eq. (22); p-weights W(M): 
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Eq. (25); raising and lowering operators J±): Eqs. (28) and (29); Littlewood-Richardson rule: Eq. (31); highest- 
weight CGCs Cm"m' : E( l- ( 36 ); normal form of highest- weight CGCs: Eq. (37); lower- weight CGCs C^'m*- Eq. (40); 
irrep index P(S): Eq. (C2); pattern index Q(M): Eq. (C7). 

To compile, type g++ clebsch.cpp -llapack -lblas on Linux, or g++ clebsch.cpp -framework vecLib on 
Mac OS X. On other operating systems, make sure that LAPACK is included in the linking process. To achieve that, 
you may have to modify the declaration of the funtions dgesvd and dgels. 



#i n c 1 u d e 
^ti nclude 
#i nclude 
^ti nclude 
^include 
^include 
^include 
#i nclude 
^ti nclude 
^ti nclude 
#i nclude 
^include 



< algor it hm > 
<casscrt > 
<cmath> 
< c st d i o > 
<cst dlib > 
<cstring> 
<f unct ional > 
<fstream > 
<iostream > 
<map> 
<numeric> 
<vector > 



// D eclaration of LAPACK subroutines 

// Make sure the data types match your version of LAPACK 

extern "C" void dgesvd. ( char const * JOBU, 

char const* JOBVT, 
int const * M, 
int const * N , 
double* A, 
int const* LDA, 
double* S, 
double* U, 
int const* LDU, 
double* VT, 
int const* LDVT, 
double* WORK, 
int const* LWORK, 
int *INFO); 

extern "C" void dgcls. (char const* TRANS, 

int const* M, 
int const * N , 
int const* NRHS, 
double* A, 
int const* LDA, 
double* B, 
int const* LDB , 
double* WORK, 
int const* LWORK, 
int *INFO); 



namespace c 1 c b s c h { 

const double EPS — lc- 



12; 



// binomial c o efficients 
class binomial.t { 

std : : vcctor<int> cache ; 

int N; 



public : 

int operator () ( int n, int k); 

} binomial ; 

// Eq. (19) and (25) 
class weight { 

std::vcctor<int> elem; 

public : 

// the N in "SU(N)" 
const int N; 



// create a non— initi ali z e d weight 
weight ( int N) ; 



// create irrep weight of given index 
// Eq. (C2) 

wcight(int N, int index); 

// assign from another instance 

clebsch : : weight &operat or — (const clcbsch : : weight &w ) ; 
// access elements of this weight (k — 1, ... , N) 
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}; 



int ^operator ( ) ( int k ) ; 

const int ^operator ( ) ( int k) const; 

/ / compare weights 
// Eq. (CI) 

bool operator < (const weight &w) const ; 
bool operator— —(const weight &w) const ; 

// element — wise sum of weights 

clebsch :: weight operator + ( const weight &w) const; 

// returns the index of this irrep weight (index — 0, 1, ■■■) 

// Eq. (Cft) 

int index () const ; 

// returns the dimension of this irrep weight 
// Eq. (22) 

long long dimension () const; 



// Eq. (20) 
class pattern { 

std::vcctor<int> elem; 

public : 

// the N in "SU(N)" 
const int N; 

// copy constructor- 
pattern (const pattern &pat ) ; 

// create p at tern of given index from irrep weight 
// Eq. (C7) 

pattern (const weight &irrcp , int index — 0); 

// access elements of this pattern (I — 1, ... , N; k — 1, ... , I) 

int ^operator ( ) ( int k, int 1); 

const int ^operator ( ) ( int k, int 1) const; 

// find succeeding / preceding pattern . return false if not possible 

// Eq. (C9) 

bool operator + + (); 

bool operator (); 

// returns the pattern index (index — 0, ... , dimension — 1) 

// Eq. (C7) 

int index () const ; 

// returns the pattern weight 
// Eq. (25) 

clebsch : : weight gct.wcight () const ; 

// returns matrix element of lowering operator J " ( I ) _ — 
// between this pattern minus M" (k , I) and this pattern 
// (I = 1, . . . , N; k — 1 , . . . , I) 
// Eq. (28) 

double lowcring_cocff(int k, int 1) const; 

// returns matrix element of raising op erator J " ( I ) _-/- 
// between this p attern plus M"(k, I ) and this p attern 
// (I = 1, ■ ■ ■ , N; k — 1 , . . . , I) 
// Eq. (29) 

double raising_coeff(int k, int 1) const; 

}; 

class decomposition { 

std : : vcctor<clcbsch : : weight> weights ; 
std : : vcctor<int> multiplicities ; 

public : 

// the N m "SU(N)" 
const int N; 

// save given irreps for later use 
const weight factorl , factor2; 

// construct the decomposition of factorl times factor2 into irreps 
// Eq. (31) 

decomposition ( const weight &factorl , const weight &factor2); 

// return the number of occurring irreps 
int s i z c ( ) const ; 
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// access the occurring irreps 
// 3 = 0, ... , size () - 1 

const clcbsch : : weight ^operator () ( int j) const; 

// return the outer multiplicity of irrep in this decomposition 
int m u 1 1 i pi i c i t y ( const weight &irrep) const; 

}; 

class indcx_adapt cr { 

std::vcctor<int> indices; 

std : : vcctor<int> multiplicities ; 

public : 

// the N %n "SU(N)" 
const int N; 

// save given irreps for later use 
const int factorl , factor2; 

// construct this index ^adapter from a given decomposition 
index.adaptcr (const clebsch : : decomposition &decomp ) ; 

// return the number of occurring irreps 
int s i z c ( ) const ; 

// access the o c cur ring irreps 
int operator () ( int j) const; 

// return the outer multiplicity of irrep in this decomposition 
int mu 1 1 i pi i c i t y ( int irrep) const; 

}; 

class coefficients { 

std:: map<st d ::vcctor<int>, double> elzx; 

// access Clebsch — Gordan co efficients in convenient manner 
void set(int factorl.state , 

int f ac t o r 2 _s t at e , 

int m u 1 1 i p 1 i c i t y _i n d e x , 

int irrcp.state , 

double value ) ; 

// internal functions , doing most of the work 
void highest_weight_normal_form(); // Eg . (37) 
void co m p u t c _h ighes t _ we ight _coef f s ( ) : // Eg. (36) 

void comput e -lower .weight _coeffs (int multip_index , int state , std : : vector <char> &donc ) ; // Eq . (40 ) 

public : 

// the N in "SU(N)" 
const int N; 

// save irreps and their dimensions for later use 
const weight factorl , factor2 , irrep; 

const int f ac t o r 1 _d i men s i o n , f ac t o r 2 _d i me n s i o n , i r r e p _d i me ns i o n ; 

// outer multiplicity of irrep in this decomposition 
const int multiplicity ; 

// construct all Clebsch — Gordan coefficients of this decomposition 

coefficients (const weight &irrcp , const weight &factorl , const weight &factor2); 

// access Clebsch — Gordan coefficients (read— only) 
// multiplicity -index — 0, ... , multiplicity — 1 
/ / factorl _ state — , ... , f actor 1 _dirnension — 1 
/ / fact or 2 -state — , ... , f actor 2 -dim en si on — 1 
// irrep -state — 0, ... , irrep -dimension 
double operator ()( int factorl.state , 

int f a c t o r 2 _s t a t e , 

int m u 1 1 i p 1 i c i t y _ i n d ex , 

int irrcp.state ) const; 

}; 

}; 

// implementation of " binomial_t" starts here 

int cleb sc h :: binomial _t :: operator ()( int n, int k) { 
if (N <= n) { 

for (cache .rcsizc((n+ 1) * (n + 2) / 2): N<=n; ++N) { 

cache [N * (N + 1) / 2] = cache [N * (N + 1) / 2 + N] = 1 ; 
for (int k = 1; k < N; ++k ) { 

cache [N * (N + 1) / 2 + k] = cache [(N - 1) * N / 2 + k] 

+ cachc[(N - 1) * N / 2 + k- 1]; 

} 
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} 

} 

return cache [n * (n + 1) / 2 + k] : 

} 

// implementation of "weight" starts here 

clcbsch :: weight :: weight ( int N) : clcm(N), N(N) {} 

clcbsch :: weight :: weight ( int N, int index) : elcm(N, 0), N(N) { 
for (int i = 0: index > && i < N ; ++i ) { 

for (int j = 1; binomial(N — i — l + j,N— i — 1) <= index; j «= 1) { 
elem [ i ] = j ; 

} 

for (int j = clem [ i ] » 1; j > 0; j »= 1) { 

if ( binomial (N - i - 1 + (elem[i] | j), N- i - 1) <= index) { 
elem [ i ] |= j ; 

} 

} 

index — — binomialfN— i — 1 + c 1 c m [ i ] + + , N— i — 1); 

} 

} 

clebsch : : weight &clebsch : : weight : : operator — (const clcbsch : : weight &w) { 
int &n — const_cast<int &>{N); 
clem — w. elem ; 
n — w.N; 
return * t h is ; 

} 

int &clcbsch :: weight:: operator ()( int k ) { 
asscrt(l <= k k <= N); 

return elem[k — 1]; 

} 

const int &clebsch:: weight:: operator ()( int k) const { 

asscrtfl <= k k <= N); 

return elem [k — 1] ; 

} 

bool clcbsch :: weight :: operator < (const weight &w) const { 
assert (w.N = N) ; 
for ( int i = 0; i < N; ++i ) { 

if (elem [ i ] — elem [N — 1] != w. elem [ i ] - w. elem [N - 1]) { 

return elem [ i ] — clcm[N — 1] < w. elem [ i ] — w.clcm[N — 1]; 

} 

} 

return false ; 

} 

bool clebsch : : weight : : operator— —(const weight &w) const { 
assert (w.N = N) ; 

for ( int i = 1 ; i < N; ++i ) { 

if (w . elem [ i ] — w . elem [ i — 1] !— elem [ i ] — elcm[i — 1]) { 
return false ; 

} 

} 

return true : 

} 

clebsch : : weight clcbsch : : weight : : operator + (const weight &w) const { 
weight result(N); 

transform (elem . beginQ , elem . en d ( ) , w.clcm . begin () , result .clem . begin () , std : : plus<int >()); 
return result ; 

} 

int clcbsch::wcight::indcx() const { 
int result — 0; 

for (int i = 0; elem [ i ] > clem [N - 1]; ++i ) { 

result +— binomial(N — i — 1 + elem [ i ] — clcm[N — 1] — 1. N — i — 1); 

} 

return result ; 
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long long clcbsch : : weight : : dimension () const { 
long long numerator — 1 , denominator — 1 ; 

for ( int i = 1 ; i < N; ++i ) { 

for (int j - 0; i + j < N; ++j ) { 

numerator *= clem [ j ] — clem [ i + j ] + i ; 
denominator *— i : 

} 

} 

return numerator / denominator ; 

} 

// implementation of "pattern" starts here 

clebsch : : pattern : : pattern (const pattern &p ) : elem(p.elem), N(p.N) {} 

clebsch : : pattern : : pattern (const weight &irrcp, int index) : 
elem ( ( irrep .N * (irrep.N+ 1)) / 2), N(irrcp.N) { 
for (int i = 1; i <= N; ++i ) { 
(*this)(i, N) = irrep(i): 

} 

for (int 1 = N - 1; 1 >= 1; 1) { 

for ( int k = 1 ; k <= 1 ; ++k ) { 

(*this)(k, 1) = (*this)(k +1.1 + 1): 

} 

} 

while (index > 0) { 

bool b = ++(*this); 



} 



asser t ( b ) : 

} 



int &clebsch :: pattern :: operator ()( int k, int 1) { 

return elem[(N * (N + 1) — 1 * (1 + 1)) / 2 + k - 1]: 

} 

const int &c 1 c b s c h : : pattern : : oper at or ( ) ( int k, int 1) const { 

return clcm[(N * (N + 1 ) - 1 * ( 1 + 1 ) ) / 2 + k - 1 ] : 

} 

bool clcbsch : : pattern : : operatorH — \-() { 
int k = 1 , 1 = 1; 

while (1 < N && (*this)(k. 1) = (*this)(k, 1 + 1)) { 
if (— k = 0) { 
k = ++1 ; 

} 

} 

if (1 = N) { 

return false ; 

} 

++(*this)(k, 1); 

while (k != 1 | | 1 != 1) { 
if (++k > 1) { 
k = 1; 
1 : 

} 

(*this)(k. 1) = (*this)(k +1,1 + 1); 

} 

return true ; 



bool clcbsch : : pattern : : operator () { 

int k = 1 , 1 = 1; 

while (1 <N&& (* this )(k, 1) = (*this)(k + 1. 1 + 1)) { 
if (— k == 0) { 
k = ++1 ; 

} 

} 

if (1 = N) { 

return false ; 

} 



24 



(*this)(k, 1); 

while (k != 1 | | 1 != 1) { 
if (++k > 1) { 
k = 1; 
1 : 

} 

(*this)(k, 1) = (*this)(k, 1 + 1): 

} 

return true : 

} 

int clcbsch :: pattern :: indcx() const { 
int result — ; 

for (pattern p(*this); p; +- |-rcsult) {} 

return result ; 

} 

clebsch : : weight clebsch : : pattern : : gct_weight () const { 
clcbsch : : weight result (N) ; 

for (int prcv = 0, 1 = 1; 1 <= N; ++1 ) { 
int now — 0; 

for ( int k = 1 ; k <= 1 ; ++k ) { 
now +— ( * t h i s ) ( k , 1): 

} 

result (1) — now — pr ev ; 
prcv — now; 

} 

return result ; 

} 

double clebsch:: pattern ::lowcring_cocff(int k, int 1) const { 
double result = 1.0; 

for (int i = 1; i <= 1 + 1 ; ++i ) { 

result *= (*this)(i, 1 + 1) - (*this)(k, 1) + k - i + 1; 

} 

for (int i = 1; i <= 1 - 1 ; ++i ) { 

result *= (*this)(i, 1 — 1) — (*this)(k, 1) + k — i; 

} 

for (int i = 1; i <= 1 ; ++i ) { 
if (i — — k) continue; 

result /= (* this ) ( i , 1) - (* this ) (k , 1) + k - i + 1: 
result /= (*this)(i, 1) - (*this)(k, 1) + k - i: 

} 

return std :: sqrt(— result ); 

} 

double clebsch:: pattern ::raising_coeff (int k, int 1) const { 
double result = 1.0; 

for (int i = 1; i <= 1 + 1 ; ++i ) { 

result *= (*this)(i, 1 + 1) - (*this)(k, 1) + k - i; 

} 

for (int i = 1; i <= 1 - 1 ; ++i ) { 

result *= (*this)(i, 1 - 1) - (*this)(k, 1) + k - i - 1; 

} 

for (int i = 1; i <= 1 ; ++i ) { 
if (i — — k) continue; 

result /= (* this )( i , 1) - (* this )(k, 1) +k~ i; 
result /= (*this)(i, 1) - (*this)(k. 1) + k - i - 1: 

} 

return std :: sqrt(— result ); 

} 

// implementation of "decomposition" starts here 

clebsch :: decomposition :: decomposition ( const weight &factorl , const weight &factor2) 



25 



} 



N( factor 1 .N) , factor 1 ( factor 1 ) , f a c t o r 2 ( f a c t o r 2 ) { 
assert ( factorl .N = factor2.N); 
std : : vcctor<clcbsch : : weight> result ; 
pattern low ( factor 1 ) , high(factorl); 
weight trial ( factor2 ) ; 
int k = 1 , 1 = N; 

do { 

while (k <= N) { 

1 : 

if (k <= 1) { 

low(k, 1 ) = std : : max(high (k + N- 1. N). high(k. 1 + 1) + trial ( 1 + 1) - trial ( 1 ) ) ; 
high(k, 1) = high(k, 1 + 1); 

if (k > 1 fofe high(k. 1) > high(k - 1, 1 - 1)) { 
high (k, 1 ) = high (k - 1 . 1 - 1) ; 

} 

if ( 1 > 1 &fe k = 1 && high (k , 1) > trial ( 1 - 1) - trial ( 1 )) { 
high(k, 1) = trial(l - 1) - trial(l); 

} 

if (low(k, 1) > high(k, 1 )) { 
break ; 

} 

trial(l + 1) += high(k, 1 + 1) - high(k, 1); 
} else { 

trial(l + 1) += high(k, 1 + 1); 
++k; 
1 = N; 

} 

} 

if (k > N) { 

result . push_back( trial ) ; 

for (int i = 1; i <= N: ++i ) { 

result. back()(i) — — result ,back()(N); 

} 

} else { 

++1 ; 

} 

while (k != 1 | | 1 != N) { 
if (1 = N) { 

1 = k - 1; 

trial (1 + 1) -= high(k, 1 + 1); 
} else if (low(k. 1) < high(k, 1)) { 
— high(k, 1); 
++trial (1 + 1) ; 
break ; 

^ trial ( 1 + 1) -= high(k, 1 + 1) - high(k, 1): 
} 

++1 ; 

} 

} while (k != 1 | | 1 != N) ; 

sort ( result . begin ( ) , result . end ( ) ) ; 

for (std::vector<clebsch::wcight>::itcrator it — result, begin (); it !— result . end ( ) ; ++i t ) { 
if (it !— r es u It . begin ( ) *it = weights . back () ) { 

-f- f-m ultiplicitics . back ( ) ; 
} else { 

weights . push_back(* it ) ; 
multiplicities .push_back(l): 

} 

} 



int c lebsch :: decomposition :: s i z c ( ) const { 
return weights . size () ; 

} 

const clcbsch :: weight feclcbsch : : decomposition : : operator () ( int j) const { 
return weights [ j ] ; 

} 

int c lebsch :: decomposition :: m u It i p 1 i c i t y ( const weight &irrcp) const { 
assert ( irrcp . N = N) ; 

std : : vcctor<clebsch : : weight >:: const-iterator it 

— std : : lowcr_bound (weights . begin () , weights . end () , irrcp ) ; 

return it !— weights . end ( ) SzSz * it = irrep ? m u 1 1 i p 1 i c i t i e s [ i t — weights . begin () ] : 0; 

} 

// implementation of " index_adapter " starts here 
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clebsch : : indcx_adapter : : index_adapter (const clebsch : : decomposition &decomp ) : 
N(decomp . N) , 

factorl (decomp . factorl . index ()) , 
factor2 (decomp . factor2 . index ()) { 
for (int i — 0, s — decomp . size () ; i < s; ++i ) { 
indices . push.back (decomp( i ) . index () ) ; 

multiplicities .push_back( decomp . multiplicity ( decomp ( i ) ) ) : 

} 

} 

int clebsch : : indcx.adapter : : size () const { 
return indices. size (); 

} 

int clebsch ::index_adapter:: operator ()( int j) const { 
return indices [ j ] ; 

} 

int clebsch:: indcx.adapter:: multiplicity (int irrcp) const { 

std::vcctor<int>::const_itcrator it — std : : lowcr_bound (indices . begin () , indices, end ( ) , irrcp); 

return it !— i n d i c c s . end ( ) && *it = irrcp ? m u 1 1 i p 1 i c i t i e s [ i t — i n d i c c s . begi n ( ) ] : 0: 

} 

// implementation of "clebsch" starts here 

void clebsch :: coefficients :: sct(int factorl.statc , 

int factor2_statc , 
int m u 1 1 i p 1 i c i t y _i n d ex , 
int irrep.state , 
double value) { 

assert (0 <— factorl.statc && factorl.statc < factorl.dimension); 
assert (0 <— factor2_state && factor2_state < factor2_dimension); 
assert (0 <— multiplicity_index && multiplicity_indcx < multiplicity); 
assert (0 <— irrcp.statc && irrep.statc < irrcp.dimcnsion); 

int coefficient-label [] — { factor 1 .state , 

factor2_state , 
m u 1 1 i p 1 i c i t y _ i n d c x , 
irrep.state }; 

elzx [std : : vcctor<int>(cocfficient_labcl , coefficient_label 

+ sizeof coefficient-label / sizeof coefficient-label [0])] — value: 

} 

void clebsch:: coefficients ::highest_weight_normal_form() { 
int hws — irrep.dimension — 1; 

// bring CGCs into reduced row echelon form 

for (int h — 0, i — 0; h < multiplicity — 1 i < factorl.dimension; i I i ) { 

for (int j — 0; h < multiplicity — 1 && j < f ac t or 2 _d i men s io n ; ++j ) { 
int kO = h; 

for (int k = h + 1; k < multiplicity: ++k ) { 

if ( fabs ( (* this ) ( i , j, k, hws)) > fabs ((* this )( i , j, kO , hws))) { 
kO = k: 

} 

} 

if ((*this)(i, j, kO , hws) < -EPS) { 

for (int i2 — i; i2 < factorl.dimension; ++ i 2 ) { 

for (int j 2 = i 2 = i ? j : 0; j2 < f act or 2 _d i me ns ion ; ++j 2 ) { 
set(i2, j2 , kO , hws, -(*this)(i2, j2, kO . hws)): 

} 

} 

} else if ((*this)(i, j, kO , hws) < EPS) { 
continue ; 

} 

if (kO != h) { 

for (int i2 — i; i2 < factorl.dimension; ++ i 2 ) { 

for (int j 2 = i 2 = i ? j : 0; j2 < f act or 2 _d i me ns ion ; ++j 2 ) { 
double x = (* this ) ( i2 , j2 , kO , hws ) ; 
set(i2, j2 , kO , hws, (*this)(i2, j2 , h, hws)); 
s e t ( i 2 , j2, h, hws, x): 

} 

} 

} 

for (int k = h + 1; k < multiplicity; ++k ) { 

for (int i2 — i; i2 < factorl.dimension; ++i2) { 

for (int j 2 = i 2 = i ? j : 0; j2 < f ac t o r 2 _d i me ns i o n ; ++j 2 ) { 

set(i2, j2 , k, hws, (*this)(i2, j2, k, hws) - (*this)(i2, j2, h, hws) 
* (*this)(i, j, k, hws) / (*this)(i, j, h, hws ) ) ; 
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} 

// next 3 lines not strictly necessary , might improve numerical stability 
for (int k = h + 1; k < multiplicity; ++k ) { 
set ( i , j , k , hws , 0.0): 

} 

++h; 



} 



/ / Gram- Schmidt orthonormalization 
for (int h = 0; h < multiplicity; +-f-h ) { 
for (int k = 0; k < h; ++k) { 
double overlap — 0.0; 

for (int i — 0; i < factor l.dimcnsion ; ++ i ) { 

for (int j — 0; j < f ac t o r 2 _d i me ns io n ; ++j ) { 

overlap +— (*this)(i, j. h, hws ) * ( * t his ) ( i 

} 

} 



j , k , hws ) ; 



} 



for (int i — 0; i < factor l_dimension ; +-f i ) { 

for (int j — 0; j < factor2_dimension; -H-j ) { 

sc t ( i , j , h , hws . (*this)(i, j. h, hws ) — overlap * (*this)(i 

} 

} 



j , k , hws ) ) ; 



double norm — 0.0; 

for (int i — 0; i < factorl.dimcnsion ; -j— (- i ) { 

for (int j — 0; j < factor2_dimension; ++j ) { 

norm +— (*this)(i, j. h, hws ) * (*this)(i, j, r. 

} 

} 

norm — std::sqrt( norm ) ; 

for (int i — 0; i < factorl.dimcnsion; -j-+i ) { 

for (int j — 0; j < factor2_dimension; ++j ) { 

j , h, hws, (*this)(i , j , h, hws) / norm ) 



hws ) ; 



} 



} 



set ( i 



} 



void clcbsch : : coefficients : : compute_highest_weight_coeffs () { 
if (multiplicity = 0) { 
return : 

} 

std : : vcctor<std : : vcctor<int> > map_coeff(factorl_dimension , 

std :: vcctor<int>(factor2_dimension , 
std : : vector <std : : vector <int> > map.states ( factorl.dimcnsion , 

std : : vect or <int>( fact or 2 .dimension 

int n.cocff — 0, n.statcs — 0; 
pattern pffactorl , 0); 

for (int i — 0; i < factor 1 .dimension ; ++i , ++P ) { 
weight pw (p.gct.wcight ()); 
pattern q(factor2 , 0); 

for (int j — 0; j < factor2_dimcnsion; -H-j , ++q ) { 
i f ( pw + q.gct_wcight() = i r r e p ) { 
map.cocff [ i ] [ j ] — n_cocff++; 

} 

} 

} 

if (n.cocff = 1) { 

for (int i — 0; i < factorl.dimcnsion; ++i ) { 
for (int j — 0; j < fact or 2 .dimension ; +-f- 
if ( map_cocff [ i ] [ j ] >= 0) { 

sct(i , j , 0, irrcp.dimension 
return ; 

} 



-i)); 



l 



j) { 

1.0); 



} 



} 



} 



double *hw_system — new double [ n _c o c ff * (factorl.dimcnsion * f ac t o r 2 _d i m en s io n ) ] ; 

assert ( hw_systcm !— NULL); 

memset ( hw_system , 0, n.cocff * ( factor 1 .dimension * f act o r 2 _d i me ns i o n ) * sizeof (double)); 
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pattern r(factorl , 0); 

for (int i — 0; i < factorl_dimcnsion; ++i , -(— (-r ) { 
pattern q(factor2 , 0); 

for (int j — 0; j < factor2_dimcnsion; ++j , ++q ) { 
if ( map.cocff [ i ] [ j ] >= ) { 

for (int 1 = 1; 1 <= N — 1 ; ++1 ) { 
for (int k= 1; k<= 1; ++k ) { 

if ((k = 1 || r(k, 1) + 1 <= r(k - 1, 1 - 1)) && r(k, 1) + 1 <= r(k, 1 + 1)) { 
++r(k, 1); 
int h — r . index () ; 
— r(k, 1); 

if f map_st ates [ h ] [ j ] < 0) { 

map_states[h][j] — n_statcs+ + ; 

} 



} 



hw.system [ n.cocff * map.states [ h ] [ j ] + map.coeff [ i ][ j 
+— r . raising_cocff(k, 1 ); 



if ((k = 1 || q(k, 1) + 1 <= q(k - 1, 1 - 1)) &fc q(k, 1 ) + 1 <= q(k, 1 + 1)) { 
++q(k, 1); 
int h — q.index(); 
— q(k, 1); 

if ( map.states [ i ][ h ] < 0) { 

map.states [ i ] [h] — n_statcs+ + ; 

} 

hw.system [ n.cocff * map_states[i][h] + map.coeff [ i ][ j ] ] 
+— q. raising_cocff(k, 1 ); 



int lwork — — 1 , info; 
double worksize ; 

double * sing val — new double [ st d : : min (n.cocff , n_statcs )] : 
assert (singval !— NULL ) ; 

double *singvcc — new double [ n _ c o c ff * n.cocff]; 
assert (singvec !— NULL ) ; 

dgesvd. ("A" , 
"N" , 

&n _c o e f f , 
&n_states , 

hw.system , 
&cn _cocf f , 

singval , 

singvec , 
&n _c o e f f , 
NULL, 

&n_states , 
&worksize , 
&lwork , 
&info ) ; 
assert f info — — ) ; 

lwork — worksize ; 

double *work — new double [ lwork ] ; 
assert (work !— NULL); 

dgesvd. ("A" , 
"N" , 

&n .coeff , 
&n_states , 

hw.system , 
&n .coeff , 

singval , 

singvec , 
&n _cocf f , 
NULL, 

&n_st ates, 
work , 
&lwork , 
&info ) ; 
assert f info — — ) ; 
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} 



for (int i — 0; i < multiplicity; ++ i ) { 

for (int j — 0; j < factorl.dimension ; ++j ) { 

for (int k — 0; k < factor2_dimension; ++k ) { 
if (map.cocff [j ] [k] >= 0) { 

double x — singvec [ n.coeff * ( n.cocff — 1 — i) + map.cocff [ j ][ k ]] : 

if ( fabs (x) > EPS) { 

set ( j , k, i. irrcp_dimcnsion — 1, x); 

} 

} 

} 

} 

} 

/ / uncomment next line to bring highest — weight co efficients into "normal form " 
/ / highest _w eight _normal_j 'o rm ( ) ; 

delete [ ] work ; 

delete [] singvec ; 

delete [ ] sing v al ; 

delete [] hw.system ; 



void clcbsch :: coefficients :: compute_lower_weight_cocffs(int multip.indcx , 

int state , 

std : : vector <char> &donc) { 

weight statcw ( pattern ( irrcp , state). gct_wcight()); 
pattern p(irrcp , 0); 

std::vcctor<int> map.parent ( irrep.dimension , — 1 ) , 

map.multi ( ir r ep _di men si o n , — 1), 

wh i c h _1 ( i r r ep _d i me ns i o n , —1); 
int n_parent — 0, n_multi — 0; 

for (int i — 0; i < irrep.dimension; ++i , -H-p ) { 
weight v(p. gct.wcight () ) ; 

if (v statcw) { 

map_multi[i] — n.multH — h; 
} else for (int 1 = 1; 1 < N: ++1 ) { 
— v(l); 
++v(l +1); 

if ( v statew) { 

map.parent [ i ] — n_parcnt++; 
which.l [ i ] — 1 ; 
if (!donc[i]) { 

computc.lowcr.wcight.cocffs ( multip.index , i . done); 

} 

break : 

} 

— v(l + 1); 
++v ( 1 ) ; 

} 

} 

double *irrep_coeffs — new double [ n.parcnt * n.multi ] ; 
assert (irrep-coeffs !— NULL ) ; 

memset ( i r r c p _c o c f f s , 0, n_parent * n_multi * sizeof (double)); 

double *prod_coeffs — new double [ n.parent * factorl.dimension * factor2_dimcnsion ] ; 
assert ( prod.coeffs !— NULL ) ; 

memset ( prod_coeffs , 0, n.parent * factorl.dimension * factor2 .dimension * sizeof (double)); 

std : : vcctor<std : : vector<int> > map.prodstat ( factorl.dimension , 

std : : vcctor<int>(factor2_dimension , — 1)); 

int n.prodstat — 0; 
pattern r(irrcp , 0); 

for (int i — 0; i < irrep.dimension; ++i , -H-r ) { 
if (map.parent [ i ] >— 0) { 

for (int k = 1, 1 = which -1 [ i ] : k<= 1: ++k) { 

if (r(k, 1) > r (k + 1 , 1 + 1) && (k = 1 || r(k. 1 ) > r (k, 1- 1))) { 
— r(k, 1); 
int h — r . index () : 
++r(k, 1); 

irrcp.coeffs [n.parcnt * map_multi[h] + map_parcnt[i]] +— r ■ lowcring_coeff(k, 1 ): 

} 

} 

pattern ql(factorl , 0); 

for (int jl — 0; jl < factorl.dimension; +- 1- j 1 , ++q 1 ) { 
pattern q2(factor2 , 0); 
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for (int j2 = 0; j2 < f ac t o r 2 _d i m c n s io n ; ++j 2 , ++q2 ) { 

if ( std : : fabs ((* this ) ( jl , j2 , multip .index, i )) > EPS) { 
for ( int k = 1 , 1 = which_l [ i ] ; k <= 1 : ++k ) { 

if (ql(k. 1) > ql(k +1, 1 + 1) && (k = 1 || ql(k. 1) > ql(k, 
— ql(k, 1); 
int h — ql . index (); 
++ql(k, 1); 

if ( map-prodstat [h] [ j2 ] < 0) { 

map_prodst at [ h ] [ j 2 ] — n.prodstatH — h; 

} 

p r o d _c ocf f s [ n_parcnt * map_prodst at [ h ] [ j 2 ] + map_parent [ : 



i - i))) { 



+= 



} 



(*this)(jl, j 2 , multip_indcx, i) * ql.lowering_coeff(k, 1); 

q2(k, 1) > q2(k, 1 - 1))) { 



1 



if (q2(k. 1) > q2(k +1, 1 + 1 ) && ( k 
— q2(k, 1); 
int h — q2.indcx(): 
++q2(k, 1); 

if ( map.prodstat [ j 1 ] [ h ] < 0) { 

map_prodstat [jl ] [h] — n _pr od s t at ++; 



prod_coeffs [ n.parent * map.prodstat [ j 1 ][ h ] + map.parent [ i ] ] +— 

(*this)(jl, j2, multip_indcx, i) * q2.1owering_coeff(k, 1); 



double worksize ; 

int lwork — — 1 , info; 

d gc 1 s _ ( "N" , 

&n.parent , 

&n_multi , 

&n_prodstat , 
irrcp_coeffs , 

&n.parent , 
prod.cocffs , 

&n_parent , 

&worksizc , 

&lwork , 

&info ) ; 
assert (info — — 0); 

lwork — worksize ; 

double *work — new double [ lwork ] ; 
assert (work !— NULL ) ; 

dgcls_ ("N" , 

&n_parent , 

&n_multi , 

&n_prodstat , 
irrcp.coeffs , 

&n_parent , 
prod_cocffs , 

&n.parent , 

work , 

&lwork , 

&info ) ; 
assert ( info = ) ; 

for (int i — 0; i < irrcp.dimcnsion ; ++i ) { 
if f map_multi [ i ] >— 0) { 

for (int j — 0; j < factor l_dimcnsion ; ++ j ) { 

for (int k — 0; k < factor2_dimension; -H-k ) { 
if (map.prodstat [j ] [k] >— 0) { 

double x — prod.cocffs [n.parent * map_prodstat[j][k] + map.multi [ i ] ] ; 

if ( fabs (x) > EPS) { 

set (j , k, multip.index , i , x); 
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done [ i ] — true ; 



} 



} 



delete [ ] work ; 

delete [] prod_coeffs ; 

delete [] irrcp.coeffs 



clebsch:: coefficients:: coefficients (const weight &irrep, const weight &f act or 1 , const weight &factor2) 
N( irrcp .N) , 
factorl ( factorl ) , 
factor2 ( factor2 ) , 
irrcp ( irrcp ) , 

factorl_dimension(factorl . dimension ()) , 
factor2_dimension ( factor2 . dimension () ) , 
irrcp.dimension ( irrep . dimension ()) , 

multiplicity ( decomposition (factorl , factor2). multiplicity (irrep)) { 
assert ( factorl . N — — irrcp . N ) ; 
assert ( f a c t o r 2 . N = irrcp . N ) ; 

compute_highest_weight_coeffs () ; 

for (int i — 0; i < multiplicity; ++ i ) { 

std : : vector <char> done(irrep_dimension , 0); 
done [ irrep.dimension — 1] — true; 

for (int j — irrep.dimension — 1; j >— 0; j) { 

if (!done[j]) { 

com p ut c_lowcr .weight _cocffs ( i , j , done); 

} 

} 

} 

} 

double clebsch :: coefficients : : operator () ( int factorl.state , 

int f act o r 2 _st at e , 

int m ult i p licit y _i ndcx , 

int irrep_state) const { 
assert (0 <— factorl_state &&; factorl_state < factorl_dimension); 
assert (0 <— factor2_state && factor2_state < factor2_dimcnsion); 
assert (0 <— multiplicity.index && multiplicity.indcx < multiplicity); 
assert (0 <— irrcp_statc && irrcp.statc < irrep.dimension ); 

int coefficient-label [] — { factor 1 _state , 

factor2_state , 
m u 1 1 i p 1 i c i t y _ i n d ex , 
irrcp .state }; 
: map<st d : : vector <int >, double>:: const_itcrator it ( 

ClZX CI ■** A f c + A . . ..^^tnv + \ / n A n ff i ^Innt 



std 



: : vector<.int >, aouDie>:: const.iterator iz { 

. find(std::vcctor<int>(coefficient_label , coefficient.label 

+ sizeof coefficient - label / sizeof coefficient_label [0]))): 



return it ! = 



} 



;lzx . end ( ) ? it — >sccond : 0.0; 



// sample driver routine 

using namespace std ; 

int main ( ) { 

while ( true ) { 

int choice , N; 

cout « " What^would ~you „ 1 i kc _ to _do?" « cndl ; 

cout « " l)_Translate^an_i— weight^S^to^its _ index _P ( S ) " « endl ; 
cout « " 2) ^Recover _an _i— -weight ^S^from^ its _ i ndcx _P ( S ) " « endl ; 
cout « " 3) ^Translate^a^pattern JVL to^its^index (M) " « cndl; 
cout « " 4) _Rccovcr_a_pattern JVLfrom _its_index _Q(M) " « cndl; 

cout « " 5) ^Calculate^C lebsch— Gordan ^coefficients _for _S^x_S ' ^— >^S ' ' " « endl ; 
cout « " 6) ..Calculate _ all _ G lebsch— Gordan ^coefficients _for _S^x_S '" « endl ; 
cout « " 0)^Quit" « cndl; 



do { 

cin » choice ; 
} while (choice < 



choice > 6); 



if (choice 
break ; 



0) { 



cout « "N„(e . g . „3) : 
cin » N; 

switch (choice) { 
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} 

case 



s 1: { 
clebsch :: weight S(N); 
cout « " Irrcp _S : ; 
for (int k = 1; k <= N ; ++k ) { 
c i n » S ( k ) : 

} 

cout « S. index () « endl ; 
break ; 



2: { 
int P; 

cout « ,; 
cin » P; 
clebsch : : 
cout « ,; 
for ( int 



Index : ; 

weight S(N, P); 

I— weight : " ; 

k = 1; k <= N; 



} 

case 



cout « ' . 

} 

cout « endl ; 
break ; 



« S(k); 



HO { 



} 

case 



3: { 

clebsch :: pattern M(N); 

for (int 1 = N; 1 >= 1 ; 1 ) { 

cout « "Row„l„=„" « 1 « " : J' ; 

for (int k= 1; k<= 1; ++k ) { 
c i n » M( k , 1 ) : 

} 

} 

cout « " Index « M. index () + 1 « endl; 
break ; 

4: { 

clebsch : : weight S(N) ; 

cout « "Irrcp^S:-"; 

for (int i = 1; i <= N ; ++i ) { 

cin » S ( i ) ; 

} 



HO { 



int Q; 

cout « "Index_(l.. dim ( S ) ) : ; 
cin » Q; 

clebsch::pattcrnM(S, Q - 1); 

for (int 1 = N; 1 >= 1 ; 1 ) { 

for (int k= 1; k<= 1; ++k ) { 
cout « M(k . 1 ) « '\t ' ; 

} 

cout « endl ; 

} 

break : 

} 

case 5: { 

clebsch :: weight S(N); 
cout « " Irrcp „S„( e . g . " ; 

for (int k = N - 1; k >= 0; k) { 

cout « ' „ ' « k ; 

} 

cout « " ) : ; 
for ( int k = 1 ; k <= N ; 

cin » S ( k ) : 
} 

clebsch :: weight Sprime(N); 
cout « " Ir r cp _S ' - ( c . g . " ; 

for (int k = N - 1; k>= 0; k) { 

cout « ' „ ' « k ; 

} 

cout « " ) : ; 

for (int k = 1; k <— N; ++k ) { 

cin » Sprime(k); 

} 

clebsch : : decomposition decomp (S, Sprimc); 

cout « " Littlcwood— Richardson_decomposition_S^\\otimcs_S ' ^— ^ \ \ o p lus _S ' ' : 
cout « " [ irrep _ index ] _S ' 5 ^ ( outer^multiplicity ) _{dimension_d_S}" « endl ; 
for (int i — 0; i < decomp. size (); +-f-i ) { 

cout « " [ " « decomp ( i ) . index ( ) « " ] ; 

for (int k = 1; k <= N; ++k ) { 
cout « decomp(i)(k) « 

} 

cout « '(' « decomp . m u It i p 1 ic i t y ( decomp ( i ) ) « " ) ^{" 
« decomp ( i ). dimension ( ) « "}" « endl;; 



« endl : 
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clcbsch : : weight Sdoublcprimc (N) ; 
for (bool b — true; b; ) { 

cout « "Irrcp^S '':-"; 

for (int k = 1; k <= N; ++k ) { 
c i n » Sdoublcprimc (k ) : 

} 

for (int i = 0; i< decomp . size ( ) ; ++i ) { 
if ( decomp ( i ) — — Sdoublepr ime ) { 
b — false ; 
break ; 

} 

} 

if (b) { 

cout « "S ' '^docs^not ^occur^in ^thc^decomposition" « c n d 1 ; 

} 

} 

int alpha ; 
while (true ) { 

cout « "Outer^multiplicity^indcx : _ " ; 
cin » alpha ; 

if (1 <— alpha && alpha <— decomp . m u 1 1 i p 1 i c i t y ( Sdoublepr ime ) ) { 
break ; 

} 

cout « " S ' ' ^docs _not _havc _ t h is _ m u 1 1 i p 1 i c i t y " « cndl ; 

} 

string filc_namc; 

cout « "Entcr^filc _namc _to_writc_to_filc_(lcavc_blank_foescrecn_output): 
cin . ignorc(1234, '\n' ); 
getline (cin , filename ) ; 

const clcbsch:: coefficients C(Sdoubleprime, S, Sprime); 
int dimS — S . dimension () , 

dimSprimc — Sprime. dimension (), 

dimSdoubleprime — Sdoubleprime. dimension (); 



ofstream os ( file_name . c_str () ) ; 
( filename . empty ( ) ? cout : os). sctf(ios 
( filename .empty () ? cout 
( file.namc . empty ( ) ? cout 
( filc_name .empty () ? cout 

for (int i — 0; i < dimSdoubleprime; — | — | — i ) { 
for (int j - 0; j < dimS ; -f+j ) { 

for (int k — 0; k < dimSprimc; ++k) { 
double x — double(C(j , k, alpha — 



os ) 

OS ) 
OS ) 



fixed ) ; 
precision (15); 

« " List ^ of ^nonzero ^CGCs^ for. 
« "Q(M)\tQ(M')\tQ(M")\tCGC v 



S^x^S '^=>. 
« cndl; 



, ^ alpha" « endl : 



if (fabs(x) > clcbsch : :EPS) { 
( file.namc . empty () ? cout 
« k + 1 « ' \ t ' « i 



: os ) « j + 1 « '\t ' 
+ 1 « '\t ' « x « cndl ; 



break : 

} 

case 6: { 

clcbsch:: weight S(N); 
cout « " Ir r cp ^S^ ( c . g . " ; 

for (int k = N - 1; k>= 0; 

cout « ' ^ ' « k : 

} 

cout « " ) : ; 

for (int k = 1; k <= N ; ++k ) { 

cin » S ( k ) : 

} 



k) { 



clcbsch :: weight Sprime(N); 
cout « " Ir r cp ~S ' - ( c . g . " ; 

for (int k = N - 1; k>=0; 

cout « ' „ ' « k ; 

} 

cout « " ) : ; 

for (int k = 1; k <— N; ++k ) { 
cin » Sprime (k); 

} 



k) { 



string filename; 

cout « "Entcr^file .nameto ..write - to - fi 1 e - ( leave - blank ^ for ^screen ^output ) : 
cin.ignorc(1234, '\n'); 
getline (cin , file.name ) ; 
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ofstrcam os ( filename . c.str ()) ; 

( file.namc . empty ( ) ? cout : os). sctf(ios :: fixed); 
( filename . empty () ? cout : os). precision (15); 

clebsch:: decomposition dccomp(S, Sprime); 
(file.name .empty() ? cout : os) « 

" Littlewood— Richardson_decomposition_S^\\otimes_S ' _= _ \ \ o p 1 us _S ' ' : 
( file.name . empty ( ) ? cout : os) « 

" [ irrep ^index ] _S ' '^(outer^multiplicity )_{dimcnsion_d_S}" « endl ; 
for (int i — 0; i < decomp . size ( ) ; ++i ) { 

(file.namc. empty () ? cout : os) « " [ " « dccomp( i ) . index () « " ] . 
for (int k = 1; k <— N; ++k ) { 

( file.name .empty () ? cout : os) « decomp( i ) (k) « ' ._. ' ; 



« endl 



lultiplicity ( decomp ( i ) ) « " ) _ { '' 



} 

for 



} 

(filename .empty() ? cout : os) « ' ( ' « decomp . 
« decomp ( i ) . dimension () « " } " « endl ; ; 

(int i — 0; i < decomp. size (); ++i ) { 

const clebsch:: coefficients C(dccomp(i),S, Sprime); 
int dimS — S . dimension () , 

dimSprimc — Sprime . dimension () , 

dimSdoublcprimc — decomp ( i ) . dimension (); 

for (int m — 0; rn < C. multiplicity; -H*n) { 

( file.namc . empty () ? cout : os) « " List ^of ^nonzero ^CGCs^ for^S^x^S ' ^— >^ 
fo 



') '); 



for (int j = 1 ; j <= N; ++j ) cout « decomp ( i ) ( j ) « ( j < N ? 
( file_name . empty ( ) ? cout : os) « " , „ alpha „" « m + 1 « endl , 
( filename . empty () ? cout : os) « " Q(M) \ tQ (M ' ) \ tQ (M ' ' ) \ tCGC" « endl; 
for (int i — 0; i < dimSdoublcprimc; ++ i ) { 
for (int j = 0; j < dimS ; ++j ) { 

for (int k — 0; k < dimSprimc ; ++k) { 
double x — double(C(j , k, m, i)); 

if (fabs(x) > clebsch :: EPS) { 

( file.namc . empty () ? cout : os) « j + 1<< ' \ t ' 

« k + 1 « ' \t ' « i + 1 « '\t' « X « endl; 

} 



break : 



} 

( file_name . empty ( ) ? cout : os) « endl; 



return 0; 
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